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1. Review of methods and results. The problem to be analysed here, one with a 
wide range of applications, is to evaluate the roots of an nth order equation 


Anz” + Aye") + +++ age’ +az+ a = 0; (1) 


the roots z may be real or complex numbers. 

There are various ways of solving the problem: 
a. In one group of methods, estimated approximations of z are substituted in Eq. (1) 
and varied by trial and error until the functional value of this equation becomes zero. 
If z is complex, the computational labor increases considerably. To alleviate it, several 
electronic mechanical computers of the analogue type have been developed, e.g. by 
H. C. Hart and J. Travis [1]. A modification of this group by W. R. Evang [2] applies 
the special concepts of the theory of servomechanisms. 
b. The well known root squaring method of Graeffe transforms the coefficients a; of 
Eq. (1) gradually into a sequence, identical with the higher squared root values. How- 
ever the serious disadvantage of this method is that it yields only the moduli of the 
complex roots, whereas the phase angles have to be determined separately by additional 
methods given by F. W. J. Olver in [3]. 
c. The process of factorization splits Eq. (1) into separate factors giving immediately 
the root values. For linear factors, i.e. real roots, the process is identical with the well- 
known Regula Falsi, Newton’s rule and Horner’s method. For quadratic factors, which 
include conjugate complex roots, the procedure is due to L. Bairstow [4]: starting with 
an estimated approximation of the root factor, the coefficients a; of Eq. (1) are trans- 
formed by long division into the coefficient set b; . The remainder terms of b; , together 
with an additional coefficient set c; derived from the b; , are used to calculate corrections 
for the root parameters. The b,- and c,-sets are computed by reduction formulae suitable 
for routine work. The method represents an extension to Newton’s rule for the complex. 
It was later systematized by W. E. Milne [5]. The procedure by S. N. Lin [6] corresponds 
analogously to a simplified correction formula; however it is convergent only for restricted 
cases. 

The above factorization process is extended here by the following results. 
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In simplifying the procedure, the transform matrix and the corrections for the root 
values are derived from the b;-sets alone. As the convergence speed depends largely on 
the first good approximations of the root parameters, a system of rules of calculated 
estimates is established. The resulting procedure is simple and converges rapidly. 

The variations of the a,-coefficients are linearly dependent on the increments of 
the root parameters and vice versa. The matrices of those nth order systems are com- 
puted by direct procedures, giving important information about the properties of the 
root systems with respect to the intended application. A special case of these linear 
systems is the error relation between the variations of the a,-coefficients and the varia- 
tions of the root parameters. 

Multiple roots represent singularities as the linear relations of the transform matrix 
cease to be valid. It is shown how the error range of the multiple roots can nevertheless 
be computed directly, making the case practically identical with that of near multiple 
roots. Numerical examples illustrate the details and the specific relations in this tran- 
sition range. 

2. Reduction formulae for factorization. A pair of conjugate complex roots p, p’, 
put into factor form as (2 + p) = 0, (2 + p’) = 0, are composed into the resultant root 
factor by multiplying 


e+pet+p)=2%+p2+D, (2) 
where’ p,; = p + p’, Po = pp’ are real numbers again. Now factor (2) is divided into 
Eq. (1) by long division 
(anz” + Aner” |) e+ aoe” + az + a) + (2? + piz + po) 

= bz"? + due”? + +++ be + (dz + Bo)/(? + piz + po). (3) 


The remainder terms b, , b) would be zero if factor (2) were already correct; otherwise 
they can be used to compute corrections for the root parameters p, , po to remove the 
remainder b, , by finally. For numerical computation the a,-coefficients shall be divided 
by a, before being processed; therefore a, = b, = 1. The reduction formulae for the 


b,-coefficient set are 


b: = —Dobise — Didier +4; , (4) 

with the remainder coefficients 
b, = pobs — pibe + a, , (5) 
bo = —Dobz +a). (5’) 


The reduction formulae b; are identical with those of W. E. Milne [5], but the remainder 
term b, in (5’) is different, as the term —p,b, cannot enter it, b, being the last term of 
the quotient. 

At first as an approximation for the general case, the correction formulae for the 
root parameters are deduced for the special case of the lowest of widely separated roots. 
According to the idea of the Graeffe method, the terms in Eq. (1) higher than 2’ can be 
neglected here: 


az’ + az+a, = 0. 
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When compared with the root factor 2’ + p,z + po , this gives a first approximation 
for the root parameters 

Po = Ao/Q2 , Di = 4,/az . 
The division form of Eq. (3) with the approximate parameters py , p, may be rearranged 
by multiplying both sides by the quadratic factor 


Ane” + +++ Ao = (dae? + +++ De? + pre + po) + be + bo = O. 


Here again the terms higher than b, are neglected according to the above idea. After 
simplifying and rearranging one obtains the new remainder 


z+ (pi + b,/b2)z + (Do + bo/b.) = 0, 
which, when compared with the starting root factor, gives the corrections 
Apo = bo/bz , Ap, = b,/be ; (6) 
they show a linear interdependence on the b, , bo . The general form would be 
Apo = koibi + koobo| 
{ ’ (7) 
Ap, = kids + kiobo 


representing a linear affine relationship with the transform matrix 


= Ko Koo 
” 7 he va (8) 
The formulae are so arranged to conform with the computation forms of Sec. 6 below, 
embodying the rules of matrix multiplication—column b,b) multiplied by row ko,koo ; 
therefore the principal and secondary diagonal of K had to be interchanged. 

The special correction formulae (6) with their matrix 

-_f[ 0 1/b 

K=[sp, = ©) 
are identical with those of S. N. Lin [6], but their application is confined to the restricted 
range of widely separated roots in consequence of the supposition made above. 

Having computed the parameters pp , p, of the quadratic factor (2), the roots can 
be resolved by the well-known formula z = —p,/2 + (pi/4 — po)'””; they may turn out 
as conjugate complex or real. In this case the two linear real factors (2 + p) could be 
alternately split off from the original equation quite as well. If the degree of the given 
equation is an odd number, this is even necessary for the one linear factor. Then the 
above general formulae have to be modified and become the well-known Horner’s 


scheme for real roots 
Division form 
Bt ++ dy = (bye + +++ be + p) + by = 0, 
Reduction formula 
b; = —pbisi + 4, (4’) 


Remainder 
bo = —pb, + dp ’ 
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‘p-correction 
Ap = kb , \7’) 
Transform coefficient 
k= 1/b:. 

Equation (7’) is plainly identical with Newton’s rule for real roots, whereas Eq. (7) 
is an analogous generalized form. 

3. Affinity transform of remainder into root correction. The remainder coefficients 
b, , bo and the corrections Ap, , Ap, are interpreted as the components of a vector in 











Fia. 1. Graphical root solution 


the b- and p-planes respectively (see Fig. 1). The problem is now to establish the trans- 
formation from one plane into the other. 

This can be done illustratively by graphical means. The starting values of b, , by 
and p; , Po are represented by the point “ or vector b“ of Fig. 1. Now two increment 
vectors Ap” and Ap” are arbitrarily chosen as probing steps, and their b-vectors are 
computed according to Eqs. (4), (5) giving the points (1) and (2). The p-increments 
may be chosen alternately along the p-axes 

(2) 


Ap” Ap;”? = 0 


Ap; = 0) Ap,” 

4 
Then (0) (1) and (0) (2) are the images of the p,- and p,-axes respectively, whereby 
the distance (0) (1) divided by Ap,” is the scale value of p, ; that of p, is obtained cor- 
respondingly. The vector (0)0 = —b” is hereafter decomposed along the p.- and 
p,-axes, and these components ratioed by the scale factors are immediately the wanted 


corrections Ap, , Ap, . Of course these computations are only applicable within the 


linear range. 
The analytical form of this approach has some advantages with regard to more 


general application. For this purpose the increment vectors are referred to the starting 
point (0) 


Ap ies p” _ p) Ab’? i pb? = aad 


Ap” - p”™ ames pf Ab‘? ii bp = | 





Or O”:_— CUD 
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These increment vectors have to comply with Eq. (7), whereby the Ab-vectors have 
to be supplemented by a negative sign, as the latter and the zero vector b® are of 
opposite sign by definition: 


Point (1) 
Ap,” = —(koAb;” + KooAbo’’) 
’ (10) 
Ap,” = — (ko, Ab,” + KooAbs”) 
Point (2) 
Api” = —(kyAb:” + kyoAbo") 
(11) 


Ap}? = —(k,, Abi” + kyoAbs”) 
Pp 


These two pairs of linear equations can be used to compute the two pairs of unknown 
coefficients ko, , koo and ky, , kio of the matrix K. According to Cramer’s rule 














| Ap” Abi” Api” Aps” 
ko = —D" = —D" 
| Apo” Ab,” Ab” Abi” ; 
Abi? Api” Ap;” Ap,” 
koo = —D" = +p" , (10’) 
| Ab\” Ap,” Ab\” Ab{” 
Ap\? Abs? Ap\” Ap\” 
ky = —D" = -D" 
- Ap,” Ab,” Ab{” Abi” ‘ 
Ab{” Ap;” Ap;” Ap\” 
kyo = —-D" = +D" ; ’ (11’) 
Ab Api” Ab ab 
where 
| Abi? Abi” Abi” Ab\” 
= | = . (12) 
| Ab\” Ab;” Abs” Ab” 





The above determinants have been transformed by mirroring at the diagonal and by 
interchanging rows. These computations are easily and rather mechanically performed 
in the special computation form of Sec. 6 below. 


The above transform matrix 
k k 
K — | 01 d 
ki kio 


for the Ap-corrections (7) is computed by three positions of the remainder b, , by . There- 
fore this approach can be regarded as a generalized adaptation of the Regula Falsi to 
complex roots. 

In order to compute the remainder vector b, , b) according to Eqs. (4), (5), the co- 
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efficients a; and b; are arranged in columns, see Fig. 2. The actual root parameters 
Po , DP; are written on a paper strip with reversed sign, multiplied and summed up with 
the other coefficients by a desk coniputer with accumulative multiplication according 
to Milne [5], as shown in Fig. 2. This procedure is used in Table 1 and works rather 
mechanically. The dotted line below 7 = 2 serves as a reminder that in Eq. (5’) for bo 
the routine term —p,b, must not be included. 























— 
Step t 2 
79 qa; 
ni aszi |b =1 \ 
: | | 
Paper Strip t -Ps b... 
(= Sa b,,, 
i a;+ t—-b; 
Ls a on, 
' a, | b, 
° a, | | b, | 


Fic. 2. bj-Paper strip plot 


4, General computation procedure with p-traces. The affinity transform of Eqs. 
(7), (8) should only be applied if the functions b, , bo of the remainder vector are suffi- 
ciently linear. If this is not yet the case, one starts with some approximate values p, , 7; , 
images the corresponding point into the 6, , bo-plane (point 1 in Fig. 3) and subsequently 
additional points 2, 3 for varying p» but p, = constant. In the }, , bp-plane they constitute 
a curve (the p,-trace), possibly showing both considerable curvature and non-uniformity 
in the 7,-scale. In the region of closest proximity to the origin 0 (point 2), one starts 








Fia. 3. p-Traces 


the next trace p,; (yp. = constant), imaging additional points 4, 5 and selecting again 
the region of closest proximity to 0. With these successive steps one can always approxi- 
mate any root value as accurately as necessary. 

During approximation one will reach a region with sufficient linearity (e.g. points 
4, 5, 6). When this happens one should switch over to the affinity transform of Sec. 3, 


as this assures the most rapid convergence. 





ws 
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As shown in Fig. 3, the images of the p-coordinate axes (p-traces) may be rotated 
with reference to the b-axes by any amount and they may have any obliqueness. But 
the sign of both systems remains the same (right hand system), following at first from 
Eq. (6) for the limit of widely separated roots, showing that both the p- and b-axes are 
hereby even coincident and, of course, orthogonal. The other limit is that of multiple 
roots, where the p-traces coincide with each other, as shown in Sec. 10 below. As all 
the actual cases lie between these two limits the above statement of the sign invariance 
of systems is proven. It can be used as a shortcut in selecting the subsequent steps 
always on the correct side of the p-traces. 

5. Shortcut procedure with calculated estimates for the p-steps. For a quick 
convergence one has to start with a good approximation of the root parameters po , 7: , 
one which already falls into the linearity range of K. A well-known approximation 
for the lowest and highest root factor exploits the idea of Graeffe’s method and neglects 
alternately the higher and the lower terms of 2‘; by taking only the first or last three 
terms of Eq. (1), one obtains 

Po = An-2 Po * ao/a, 
G-Max. asa, = l, G-Min. . (13) 
Pi = Gn-1 y= a,/d2 
These formulae are fairly accurate if the root values are sufficiently separated. If however 
they reverse their turn or contradict each other, the roots must lie close together. Another 
approximation is furnished by the mean root PomPim - By multiplying the h = n/2 
quadratic factors one obtains 


(2° + put + pore" + Pre + Por) +++ @ + Piz + Por) = 0, 
or 27+ } 2 Pu” + *** PoiPo2 *** Pon = O. 


Comparing with Eq. (1) 


> pu = d,-1; TI pon = 4), 
Dom -_ (T] Por)” eee a. 


Pin = 1/h DY) pur = Gy-1/h, 


The normal procedure is always to eliminate the lowest factor in the successive elimination 
process. If the G-values of Eq. (13) get entangled, the minimum values can be derived 
accurately enough from the mean value and the maximum value ‘already known by 
G-Max. for h = 1. These shall be termed as Estimated Min. (£-Min.) to distinguish 
them from the Graeffe Min. (G-Min.). 

E-Min. 


and by definition 


(14) 


Poxia = Dom /Poress = Pom — am seni Pom)» | 

(15) 

Pimia ~ Pim — (Distas ite Pim): 

In numerical computation, the values (13), (14), (15) are compiled in tables (see Esti- 

mation Section in Table 1). The H-Min.-values are then taken as the starting steps in 

the computation section of the tables; these values are already rather good approxi- 
mations as proven by experience. 
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The next step 2 is taken on a po-trace with a roughly estimated Ap, ; the b-remainder 
vector of step 1 is plotted in a b, , bo-diagram together with a fictitious step o with 
Po = 0, p; = 0, which gives b, = a, , bo = ap ; the origin is projected on (0 1) as 0’, 
and by linear interpolation one takes roughly App = po(1 0’)/(o 1), including the correct 
sign (see Fig. 4). After having the remainder vector of step 2 subsequently transferred 
to the diagram, one can already roughly plot a first pair of pop,-traces according to Sec. 4, 
from which the corrections for step 3 can be derived (Fig. 4). After having thus computed 


tb, 











Fic. 4. Remainder-root-diagram 


3 points as preliminary steps near the origin and obviously within the linearity range, 
the transform matrix K can now be determined (see Table 2). 

6. Numerical example, convergence number. The above outlined method will be 
further clarified by means of a numerical example. For convenience of comparison the 
algebraic equation to be solved is taken from Lin’s paper [6] (Hitchcock’s example) 


2° — 3.0122’ + 3.2252° + 1.021z° + 6.9862* — 21.8872° 
+ 8.1102? + 5.901z + 23.889 = 0. 


It is assumed that the a,-coefficients are correct to 1/2 unit of the last given digit. 
Table 1 gives the solution of the lowest root couple (h = 1). It consists of the esti- 





Estim. | G-Max. G-Min.|Mean E-Min. 
Po 43.2 +3 | 2.2 | +! 
P, -3 +1 -1 





+! 














Preliminary Steps ——e=— K - Steps —+ 








Compu-S*F|_! 4 8 i es 1 8 
tation Pp.l+i 70.9 = |+0.9 |+0.861 |+0.8605 

ato +1.0 0.9 § |+0.9% 0.9430 
i a, Pp b; \* i+ 0.945 |+0.9 





8 oe +1 +1 lea let +1 
T\-3.012 |-4%.012 [- 4012 |-3.912 |~3.957 -3.9855 
6 |+3.225 + 6.237 |+6.337 |+s.846 |+6.103 |+6.094 
§ |#1.021 [- 1.204 [-1.705 [-0,.720 |-1.339 |-1.322 
& |+6.986 |+ 1.953 + 2.988 |+2.373 |+2.997 |+2.989 
3 21.887 |~22.636 |-23.341 |-23.37S |-23.566 |-23.568 
2 |+ 8.110 [428.793 [428.762 |+27.012 |+27.799 |+27.763 
1 [45.901 |- o.28e |- 1.856 |+ 2.628 |- 0.079 |+0.001 | 
© }+23.889 |- 4.90% |- 1.997 |-0.422 |-0.04%6 |-0.00) 
f (0.6) (0.9) 0.03 0.02 - 





























TABLE 1 Rooth = 1 
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mation- and computation section. After the preliminary steps have been set up, the 
K-matrix is computed in Table 2 according to Eqs. (10’), (11’), (12); the sequence of 
operations is run mechanically along its symbolical sequence-circuit. By some K-steps 
more the remainder vector b, , bo is finally put to zero or below the error limit. The 
process converges rapidly. 





' 
Sequence- 
Circuits, Step ' ¢2) 3 
<< 


== 4p, | +°.1 oO 
r—tszap,| o | -o0. 
ab, |+1.see8 [+4482 

| 4b, |-2.907 |+1.575 


——D,* | + 18.55 
K Ke, | Keq|-0.0101 [+ 0.0288 
k bk,» 0.0187 |+0.0103 


TABLE 2 Transform matriz K, 





















































After the lowest root is thus eliminated, the b,-coefficients of the last step are taken 
as the a,-coefficients of a new equation whose degree is reduced by two. These new 
abridged coefficients are again processed according to the outlined method, eliminating 
the next lowest root. This process is continued, until the last root factor has been ex- 
tracted. 

Obviously the abridged a,-coefficients already contain some errors from the former 
root factors, which reflect again in additional errors of the new root parameters. In 
fact this error propagation is rather small and may usually be disregarded. But for a 









































h ' 2 3 & 
i+ a; bin 
error aa t coos 
8 i+! +1 +1 +! +1 
7/-3.012 |-3.9ss -5.098 +0,007 +0. 010 
6/+3.225 [46.09% 411.623 +0.591 -L426 
5 |+1.021 -b322 -tt, 831 +2.786 -3.335 
& |+6986 [+2.989 + $.697 +13.827 +3.583 
3 21.887 |-23.508 — 7.3% + 12.45% +%. S53 
2|+8.10 |+#27.763 +10.69% + 8.993 + 5.096 
- te — — + — — errorabli— — — — — — — — ee — — — -4 
1 |+S5.901 [+0001 *.001 |-o0.00% *.008 |-9.018 + .°19 |-o.012 *.058 
O 23.889 |-0,001 2.00) |40002 *.0:3 [40.013 *7.0%2 |40,034% =.!e% 
error ap - 
Ppp] * 0.8605 *.c0003/+ 2.2337 *.0001 |+ 2.655 = 0006/4468) +.0009 
P,p|+ 0.9430 *-00002/+2 0g55 %.0001 |-3 019 %.000%/-3 022 +t.0002 
K kKo,| Koo} —.0101 |+.0288 | -.0017)|-.0102 |+.0%41 |+.03827 | -.0177 |-.0010 
h ky, | Kio] 4.0187 [4.0103 | —.0114 |+.000% |-.0176 1-.0168 |+.010414+.0040 








TABLE 3 Unabridged b;-coefficients, p-root parameters and K-matrices 


closer analysis this influence can be eliminated by going back to the original a,-coefficients 
for each root factor in a second round, computing the unabridged 6,-sets, K-matrices 
and root corrections Ap as above. These results are compiled in Table 3 and moreover 
are later required for the system analysis of Sec. 8 below. 
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The convergence of the computation can be checked by the “convergence number’’ 
defined for the gth step as 


B. = | bifo” |/| bi% |, (16) 


i.e. the ratio of the absolute amounts of two subsequent remainder vectors. For the 
first K-step the denominator is the mean of the | b,,, | of all three preliminary steps as 
they all contribute equally to it. The convergence number is characteristic for the 
deviation from linearity, the error AK of K and the error A(Ap) of the correction Ap, 
namely 

B = AK/K = A(Ap)/Ap. (16’) 


For the preliminary steps the following order of magnitude should normally be attained 


mee coed, 
and for the K-steps 


Bx = 0.01 --- 0.1. 


The convergence of the K-steps is very rapid; for example 8, = 0.01 indicates a gain of 
two digits per step. It is remarkable that these limits are usually reached by means of 
the first K-matrix, as for example in Table 1. If this should not happen, however, one 
has.simply to compute the subsequent K-matrix. 

It is instructive to compare the exact results of Table 3 with those of Lin [6], according 
to the simplified form (6’); only the lowest and most separated root couple h = 1 has a 
sufficient overall accordance between the simplified K and the complete form to give 
convergence. For the other roots the secondary diagonals of K take over in importance, 
or even the signs of the principal diagonals reverse, thus changing bluntly from conver- 
gence to divergence. In order to overcome these difficulties, Lin has to increase the 
separation of the different roots by supplementary methods, which makes the procedure 
tedious and the results inaccurate. This lack of convergence is also investigated in a 
paper by F. B. Hildebrand [7]. 

7. Synthetic multiplication. As a final check all the extracted root factors are 
multiplied by each other; the resulting coefficients should be identical with the original a; . 
The computation is compiled in Table 4 and is executed by means of a paper strip plot 
according to Fig. 5. As a rough estimate the mean relative deviation of the a; may be 
































hi » | a | 2 1 | 

Po} +4681 |+2.6S55 | +2.2337/+0. 8605 original aa, /o, 

P,|-3.022 |-3.019 |+2.0855/+0.9430 Coeff. a; 1o* 
i 86 i +1 m +1 - 

7 -3.012 |-3.012 ° 

6 i+ +3.225 |+3.2285 ° 

s 13,955 |+1.01Q9 |+1.021 | 20 

by los |+6.09% |+6.983 |+6.986 | 4 

3 6.041 \-1. 32m -21.891 |-21.887 2 

2)! +16.489 lez 988 |+ 8.106 |+8.110 s 

1 | -3.022 |-22.18s \waoun + 5.897 |+5.901 7 

O | +4 681 lora.ves |*27. 760 23.887 |+23.889 ' 











TABLE 4 Synthetic multiplication 











e 
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taken as the mean relative error of the root parameters, namely 

Ap/p = [Aa;/a;Jmean = 5° 10™*. 
This already is in good accordance with the results of the exact error analysis given in 
Sec. 9 below. 





“? b; oes Q= 


n 














3 ‘ 1 
Paper a +P=- b; Qa; 
Strip + b 

‘ ot os ims @;., 


+1] b, 











Fia. 5. a;-Paper strip plot for synthetic multiplication 


8. Variations of a,-coefficients and p,,-root parameters. Variations in the a,-co- 
efficients of the algebraic equation naturally produce corresponding variations in the 
Pis-root parameters and vice versa. The first letter d of the double subscript of pz, 
stands for the degree in the quadratic factor and is alternately 0 and 1, the second h 
stands for the number of the factor. For small variations there is a linear relationship 


th 
Aa; = Z. Q; anAPan ’ (17) 
dh=01 
or by inversion 
n-1 
Apan = = Pan ;Aa; ° (18) 
i=0 


The coefficients a;,, and p,,; are the matrix elements of the linear systems, and the 
objective is to compute them numerically. Those matrices are very useful for application, 
as they give the interdependence between changes of the root parameters or motion 
characteristics p,, and the physical constants a; of a dynamic system for example. 

To deduce the matrix a,,, of Eq. (17), the original Eq. (1) is set up with the Ath 


quadratic factor 


D az’ = Do dye’ 72? + pirsz + Por). (19) 
7=2 


1=0 
The indices i and j have here to be discriminated because they start from different 
origins. By differentiation of a; and p., term by term 
Maz’ = bjz'~*Apos , 
after shortening with 2’ = 2’ and comparing with Eq. (17) 
Gin = On, *=j-2=0---n— 2. (20) 
Differentiating a; and p,, term by term 
Maz’ = bz’ "Ap, , 
after shortening with 2’ = z’* and comparing with Eq. (17) 
Qi. = ba, t#=j-l=1---n-—1. (20’) 
Thus the matrix elements a,,, are identical with the b,,-coefficients of Table 3, only 
lowered in their degree by 2 and 1 respectively. For completeness the a,4,-matrix is 
compiled in Table 5. 
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The matrix elements p,, ; of Eq. (18) can now be derived by simply inversing the 
matrix of Eq. (17) according to the well-known algorithm of Gauss-Cholesky-Bana- 
chiewicz [5] and [8]. This has actually been done for this example with sufficient accuracy. 























dh| o1 | ut 02 12 | 03 | 13 O% 1% 
i Fae ET Py Ta re lee ° +1 
© |+! ee +! |-5.098 }+1 |+0.007 | +1 + 0,010 
5 |-3.955 |+6.00% |-5.098 |+11.623 |+ 0.007 +0591 |+0.010 |-1.426 
& 146.094 a 322 |+11.623 \-u.e34 + 0.591 |+ 2786 |- 1.426 -3.335 
3 |-1.322 |+2.088 -0.831 |+5.607 |+ 2.786 [3.827 |-3.335 |+3.s03 
2 |+2.989 [23.se8 |+ 5.607 |~7. 341 lemear [riz we +3.583 [+4553 
1 23.568 27.763 [-7.34) [410.694 [+i2.454 + 8.993 |eu.553 [+5096 
O 27.763 ° +10.69% | © | 8.993 | ° +5.096 Y 








TABLE 5 Matriz a; ar 


But the labor for solving this 8th order system is already rather high; it could be alle- 
viated eventually by using automatic calculating machines, suitable for this type of 
problem, like that by R. R. M. Mallock [9]. 

However, there is another and direct approach to this problem working even more 
quickly and accurately. The coefficient a; of Eq. (1) is given an increment Aa, = 1. 
At first its effect b,, ; , bo, ; on the remainder vector for the specific root factor h is com- 
puted according to the reduction formulae (4), (5) and the paper strip plot of Fig. 2; 
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k ° |b, = ~ Poi... -P, b,.., 
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° o | b= -p, 6, = b,, ; 








TABLE 6 Singular formation of increment vector ban ¢ 


this leads to the formation of the increment vector b,, ; for the order 7 (Table 6). Now 
the values b,, ; are built up step by step beginning with 7 = 0, and one sees that they 
can be telescoped into one compound for 7 = 0 ton — 1, as in Table 7; there is no incre- 














TABLE 7 Telescoped array of increment vectors ban « 








i bi, i | beh; 

) re) +1 

i +1 ° 

‘ Ph bi, ~Poh by, ' 

3 Poh Ons ~Pih ihe |~Poh ihe 

i Por ©ihi-2 “Pin Oh i- |~Por Outn i-t 


























1955) ROOT SYSTEMS OF ALGEBRAIC EQUATIONS 13 


ment for a, , because it is assumed that a, = 1 and constant. The computation is: con- 
veniently made along a paper strip plot, as in Fig. 6. Thereafter the increment vector 


i-2 i- Paper Strip 




















by; ° +1 -P, 


ohi 





























Fic. 6. Paper strip plot for bg, ;-increment vector 


bz, ; has simply to be transformed into the increment vector p., ; by means of the transform 
matrix K, already compiled in Table 3 for the original a; 


Pari = K + bans ’ (21) 


symbolizing the matrix multiplication of Eq. (7). The numerical computations are 
compiled in Table 8. 
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h Poh Dani , 2bi,; Banm 
Pih Pdhi| 0° ! 2 3 ho s 6 7 Z pha: Pdhm 
i |togeos |d to [+t aqua booze +o78s Lo76% to.04s tools | 3.488 | 1.87 
|+ 0.9430 o| +1 © |-0,861 |+0.8N |-0.024 +0.675 0.657 -0.039 | 3.288 | 1.81 
—,0101 |+,0288 © |+0.029 0.010 |-0.016 }+0.023 }+0.008 |-0.012 }+0,019 |-0.007 |.0023%/ 0.0%9 
+ .0(67|+.0103 | +0.010 # 0.019 |-0.027 +0.009 +0.01% |}+0.021 #0008 40.011 |.00209| 0.04%6 
2 |+ 2.23% 1] Oo +1 -2.086 |+2.117 |+0.24% |-5.238 |+10.38 |-9.955] 244.2] 15.6 
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|- 3.019 o} +t © |-2.085 |~8.015 |-17.15 |-30.49 |-46.53 |-S9.s0] 7002 | 83.7 
+.0%%1 |+.0327 © lro.031 +0,042 |+0,0%% |+0,022 |-0.052 |- 0.21% |-0.508 |-0967 | 1.247 | 1.12 
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TABLE 8 Increment vectors bani , Pans , error factors ban m , Dahm 


For completeness the above deductions have to be modified also for linear factors 
(z + p). Therefore Eqs. (17), (18), (19) have to be supplemented respectively by the 
terms 


Deandp, (17), Depsda; (18), De bye’ "@ +p). (19’) 
by differentiating Eqs. (19), (19’) term by term with respect to a; and p, 
Aa,z* = bjz’""Ap, , 
after shortening with z' = z2’* and comparing with Eqs. (17), (17’) 


Qn = ba, t=j-—-1=0---n-1. (20’) 
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The reduction formula for the b-increment number b,,; of the remainder b, is, in analogy 
to the formula of Table 7, 


ba: = — Dron i-r 5 with bo,o = +1, t~=0---n-1, (22) 
from which the p-increment number p,; is derived by 
Pu = kon « ’ 


where k is the known transform coefficient of the linear root h. 

9. Error analysis. The errors of the a,-coefficients are Aa; = +e. The objective 
is to deduce the errors Ap,, of the root parameters according to Eq. (18). With regard 
to the error propagation this means that the terms p,,; have to be summed up for current 7 
as an RMS-value. The errors Ab,, of the b-remainders are formed analogously 


n-1 1/2 
Abs, = tebyr m . Din am | b> vs, | ) ) 


f with the b-error factors f (23) 
Abon - sbebon m ’ Don a= [ Zz bon * hae 


and 
ADon = +CDorm J Porn = [ b> wT aid 


with the p-error factors . ( . 
ADin = LePinm » Pim = [ pO ae 


In the case that the errors e are not equal, it would be quite as easy to sum up the terms 
of Eqs. (23), (24) with the specific and respective weights. The data required for the 
numerical computation of the error factors bz, » and pz, » are already contained in Table 8; 
only the columns for those RMS-values have to be completed. Now Table 3 can be supple- 
mented with the numerical values of the b- and p-errors. It is interesting to see how, 
in the course of the successive steps, the remainder vector b,, finally falls within the 
range of the b,,-error. From then on, no additional correction Ap,, would make sense. 
It can be seen from Table 3 that the b,,-remainders of the last steps are lying just inside 
the b,,-error range. The p,,-errors as established according to Eq. (24) show that the 
calculated root parameters are correct to the last given digit. 

10. Multiple and near multiple roots. Multiple roots offer interesting singularities 
which are best understood if one studies the transition of near multiple roots to the 
limit of multiple roots. 

After having separated from a given algebraic equation all root factors by the above 
normal procedures—starting from the lowest and highest root alternately—at last a 
residual equation remains whose root factors may lie close together. The following 
quartic equation is given as a numerical example to demonstrate the approach 


z* + 4.3162" + 10.0352” + 11.6052 + 7.230 = 0. 


(24) 


In a first step the mean root M according to Sec. 5, Eq. (14) po = ao” = +2.689, 


Pp, = a;/2 = +2.158, is taken as a probing factor. The remainder vector b, , bo is com- 
puted according to the procedure outlined above and results in b, = —0.000802, 
bo = —0.000818. It is already smaller than the last digit of the given a,-coefficients so 


that the b,-computations have to be made to the 6th digit. If the b, , bo-remainder 
lies within the b-error area, the mean root can no longer be dissolved into separate 
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factors, thus constituting physically and practically a double root. To decide this question 
the remainder vector has been plotted in the b-chart of Fig. 7 together with the b-error 
area pertaining to an error e = +0.0005 of the a,-coefficients according to Sec. 9; it 


follows that 
buen = 3.09, Dos = 6.47, 


002-1 











° 
ff 
Oo 
iP. * Oot .002 ood +b, 
" fz 
tp, = pom) 


+p(O,) 








+P(O,) +R,(0,) 


Fia. 7. b-Chart 


ptm) 


and the eb,,,, are taken as the principal axes of the error ellipse. One sees from Fig. 7 
that the remainder vector b, , bo of point M (step 1) actually lies already within the 
error area. 

For further clarification the details of the b- and p-charts have to be imaged into 
each other reciprocally. First a series of po , p,-parameters are chosen along a set of 
Po- and p,-traces computing turn by turn the b, , bo-vectors by means of the b,-compu- 
tation tables. This furnishes the images of the p-traces in the b-plane; their intersection 
with the b-error ellipse is interpolated back into the p-plane, giving there the actual 















‘ul +P(m) P,(0;) 
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we we es 
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i; oe P- error erea 
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Fia. 8. p-Chart 
p-error area; it is oddly shaped, but with central symmetry around the mean root point 
M, which obviously proves that M stands more generally for multiple roots also; it is 
an area overlapped several times, while the outer curve envelopes the different parts. 
Characteristically the p-error area is considerably enlarged compared with the b-error 
area by a factor of approximately 20 for this example. 
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The images of the p-traces in the b-plane show most impressively the singularity 
of point M (Fig. 7). The negative and positive portion of the p.(M)- and p,(M)-axes 
respectively fold together into straight lines touching themselves only in M with an 
infinitely large curvature instead of intersecting each other as in the normal case. Thus 
the component rays of Fig. 1 for example become indefinite, and all the procedures 
based on linear interdependence become inapplicable for this point M. Traces with 
central symmetry in the p-plane, e.g. the po. , p.-axes through 0, , 0, respectively in 
Fig. 8, have identical images in the b-plane, but with reversed sign. The reason for this 
singularity is that in the relationship between the b- and p-vector in the neighborhood 
of M the linear terms are missing, leaving a quadratic form in pp , p, 

b, =¢ + C1000 + Ci .o1PoP1 + C1191 7] 
(25) 
bb =o + Co. 00Po + Co.o1PoP1 + Co.11P ) 


b, , bo and po , p; denoting vector components referred to point M. If the signs of a pair 
of parameters pp , p; are reversed, indicating a central symmetrical p-vector, the values 
b, , bo (b-vector) remain unchanged, in accordance with the above Eq. (25). Since the 
linear terms are missing, consequently the determinant D in Eq. (12) and b in Eq. (7) 
tend towards zero rendering the transform matrix K indefinite and Ap in Eq. (7) in- 
determinate. 

The next problem is to dissolve the mean root into its separate factors. But it can 
be solved only if the a,-coefficients are known more accurately; it shall be assumed 
now that they are accurate to the 6th digit beyond the decimal point (e = +0.0000005). 
Again the b;-computation procedure is applied. But the shortcut of Sec. 5 has now to be 
modified for the transition from the singular point M to the linear range of the root 
point 0 (Fig. 7). As stated above, the linear terms in the neighborhood of M have vanished. 
Eq. (25) is simplified for a first rough estimate 

by = cpu or Pu = Coz’, (26) 


where by , py are the small vectors from M in the b- and p-plane respectively. Considering 
transitorily for the lower root h = 1 also a step o with pp = p, = 0, bd; = a, , bo = a 
and interpolating the vector b(Mo) linearly by the perpendicular 0’ from 0 (see Fig. 7), 
one gets a first correction p(M2) for the lower root by applying Eq. (26) 
p(M2)/p(Mo) = [b(M0’)/b(Mo)]'”’. (26’) 


As the next step 2 is planned along a p,-trace, one takes the po-components instead of 
the p-vectors as a rough estimate; numerically 

p(Mo) = —p.” = —2.689, p(M2) = Ap.” = —0.03, 
where the superscripts (1) and (2) stand for step 1 and 2 respectively (see Fig. 7). The 
next steps 3 and 4 are taken in a differential form along a p,- and a pp,-trace respectively, 
whereby their directions are estimated according to Sec. 4 and the step increments are 
derived by logarithmic differentiation of Eq. (26) 


Ap/pu = Ab/2by e (27) 


A view of Fig. 7 and a check of the convergence number 8 show that the preliminary 
steps move steadily towards the root point 0. Now the transform matrix K, can be 
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computed from the last three preliminary steps (h = 1). As a consequence of the excessive 
curvature around M, the convergence rate is rather slow, and one has to improve K 
for two steps more. The computation is concluded with step 8, where the b-remainder 
drops beyond its error range. 

The parameters of the higher root h = 2 are identical with the b,-coefficients of the: 
last step for h = 1. But as a check these root parameters can also be set into their owm 
b,-computation table, and it follows that the b-remainder is again zero within the cal-- 
culated error limits. The resulting root parameters and K-matrices are compiled in: 
Table 9. Finally those matrices are used to compute the pq, ,-increment vectors and the: 
Pan w-error factors. The errors Ap for e = +0.0000005 are included in Table 9. 





h ' 2 


vt” error ap 
10°" Poh [t2.04757 *-0000% |+2.7306/ *.000 0% 





Dim 2.13957 7.00002 42.17643 *.00002 


k,, | K,,|+20.6)+05 26.7 |-0.1 
K K,, [K,, |+22.1 | -48 22.1 |+9.6 




















TABLE 9 Multiple root characteristics pa, and Ky 


For the transform matrix K there exists the following remarkable relation within 
the computational error limits, as verified by Table 9 


K, = —K, ; (28) 
meaning that all elements of the matrix reverse their sign. This theorem can be exactly 
proven by differentiating Eq. (25) 


Ab; = (c; .002Po + C1.o1P:1) Apo + (C1.01Po + C1.112p:) Ap, ’ 


Ab, = (Co. 002Po + Co.o1P:1) Apo + (Co.01Po + Co.112p:) Ap, : 


The parameters ~o, , Pi, for the root values h = 1 and 2 (referred to M) differ only in 
SIZN Poo = —Poi » Pi2 = —Pi - Consequently the matrix elements of the above equation. 
set also reverse sign for the root values h = 1 and 2, and this property still holds true 
for the inversion matrix K, q.e.d. The analogous relation obviously exists also for multiple: 
linear factors between their transform coefficients k 

ke = —k, . (28’) 
The above K-relations can be used either as computation check or as a shortcut for the: 
computation procedure. 

At last Fig. 8 is supplemented by the p-error ellipses for the original e = +0.0005. 
of the a,-coefficients, in order to compare them with the actual p-error area. Although 
the linear range is considerably exceeded thereby, the linear error areas still roughly 
cover the actual area; they will naturally flow into each other as the matrix K tends 
at M to infinity. 

The transition of near multiple into multiple roots can now be clearly seen in Fig. 8; 
the discreet roots 0, , 0, move simply towards M deforming gradually their error area, 
until it reaches finally a double symmetric shape with reference to the p.(M)- and 
p,(M)-axes. Physically and practically there is no difference between multiple and 
near multiple roots, as their large p-error area covers their small variations. 
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SEICHE IN RECTANGULAR PORTS* 


BY 
JULIEN KRAVTCHENKO 
(Université de Grenoble) 
AND 
JOHN 8S. McNOWN 
(State University of Towa) 


Introduction. In an extension of a recent publication by one of the writers [1],f 
the excitation by an external wave of the water within a harbor is studied for the case of a 
harbor which is rectangular in plan and for which the walls are vertical and totally 
reflecting and the depth is constant. The pattern of seiche, or mass oscillation, within 
the rectangular basin is predictable for a given position and size of the entrance and 
for known characteristics of the external wave. The result is useful in the interpretation 
of occurrences both in actual harbors and in laboratory models. 

The goal of the analysis is the derivation in series form of a velocity potential from 
which the periodic velocities and surface displacements can be calculated, the period 
of the motion being that of the exciting wave. A special boundary condition for the 
entrance must be devised from observations in the laboratory. A solution to the Laplace 
equation satisfying the mixed boundary condition is obtained from an expansion in 
terms of eigen-functions and auxiliary functions, as outlined previously in note form 
[2]. Solutions for two distinct cases must be obtained depending upon whether or not 
the imposed period coincides with that for one of the eigen-functions; these are designated 
as resonant and non-resonant, respectively. Similar results were found in [1] for a harbor 
of circular shape. From that study, a few elements essential to the following develop- 
ment are repeated herein. The authors wish to acknowledge the useful suggestions 
during the preparation of the manuscript made by Professors C. 8. Yih and E. N. Oberg 
of the State University of Iowa, Professor J. Kuntzman of the University of Grenoble, 
and Mr. A. Apté of the Poona Laboratory in India, at present a graduate student at 
the University of Grenoble. 

Statement of the problem. Cartesian coordinates are utilized, with the origin in the 
undisturbed free surface and the z-axis directed vertically upward. The depth of water 
is h and the dimensions of the harbor are a and b, so that the undisturbed liquid occupies 
the space — h < z < 0;0 < x < a; and 0 < y < Bb. The opening or entrance to the 
harbor, which extends over the entire depth h in one of the walls of length b, is defined 
by ~-h<2z;0<a<2< 8 < a;y = 0. The free surface is denoted for convenience as 
D and its periphery as I. 

The periodic wave causing the mass oscillation within the harbor is taken to be a 
plane Stokian wave of period 7’; it approaches the harbor over a horizontal bottom which 
is at the same elevation ( = — h) as that of the harbor. The mechanism whereby the 
excitation within the harbor is caused was described in the first paragraph of [1]. To 
complete the statement of the problem, the velocity component V,(z, z, ¢) normal to 
the plane of the entrance must be defined, ¢ denoting time as usual. The variation of V, 
can be determined only from experiment. 


*Received Feb. 20, 1954. 
tReferences indicated by numbers in brackets are collected at the end of the paper. 
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Fortunately, for the occurrence being studied, the variation of V, can be closely 
approximated from the assumption that a clapotis, or standing wave, forms in the 
immediate vicinity of the entrance (see Chap. III of [1]), the amplitude being twice 
that of the generating wave. Hence, from the classical theory of Stokes, 


AV, re 2nt cosh [k(z + h)] , 
; =F aa” (1) 








V, = 


where k is the reciprocal of a reference length, V, is the amplitude of the velocity varia- 
tion, and A is a constant which depends upon the phase of the clapotis at the entrance. 
It is evident that the velocity is discontinuous at either end of the entrance. 

Within the harbor the motion is restricted to sinusoidal oscillations of the same 
period as the external wave and ‘is assumed to be irrotational. The liquid is, of course, 
taken to be incompressible. Definition of a velocity potential y(z, y, z, t) is the goal of 
the analysis. As shown in [1], 

e = (V,/A)F(z, y) , (2) 


where F(z, y) is an analytic function appropriate for the harbor plan in question. As 


V’o = 0, 
V’F+kF =0. (3) 


The expression for g has been chosen so that 


dy | 
nit = 0 
on = 


is satisfied; the condition of Poisson for the free surface can be expressed in the form, 
4x’ /T” = gk tanh (kh). (4) 


The boundary conditions must include the statement that the velocity normal to 
the solid walls is everywhere zero, i.e., dF /dn = 0 along the intersection of the plane of 
the free surface and the solid walls. The only boundary condition remaining is that for 
¢ in the plane of the entrance. 

To be rigorous, continuity of the two velocity fields immediately inside and im- 
mediately outside the entrance would be required, although this stringent requirement 
is not entirely essential from the physical point of view. Hence only the necessary con- 
dition that the velocity potential give the same normal velocity as defined by (1) is 
imposed. From (2), this condition becomes . 

dF/dn = A, (y¥=0,a< 2 < 8), (5) 
where n is the inward-drawn normal. The problem is thus the definition of F in such a 
manner as to satisfy (3) and the conditions imposed. 

It is essential to demonstrate that the boundary-value problem is properly posed 
and to determine whether a solution is possible. To this end, it is noted that the value of 
k is known from (4) if values are assigned to the depth h and the period 7 of the oncoming 
wave. The quantities k,,,, and ®,,,, (n, m = 0, 1, ---, ~), which are, respectively, the 
eigen-values and eigen-functions of (3) for the free surface with the boundary condition 
dF/dn = 0 around the edge I, are essential to the development. It follows that 


T"@.0 + a - 0, (6) 
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and d®,,.,,/dn = 0 along I’. These values and functions are as follows: 


2 2 
ee | ae. 
 - = ¢ (" + ", 
= __2 ner mary 
Po 0 —_ (ab)'”” ’ #,,.m — (ab)'” cos a cos b 5] (7) 


“ nex (2)" mry 
2.0 = (5) s * Po.m = oe |— 


If do is an element of D, 


/ l i d= a, (8) 


in which the double subscript n, m has been replaced for convenience by i, and 6;; is 
the Kronecker delta. The uniqueness of this result is established in classic reference 
works such as [3]. 

Were it not for the singularities at either edge of the entrance, the existence of a 
solution could be readily demonstrated. For a simple approximation to reality, however, 
V,, was assumed to be independent of y (and hence discontinuous), and for this case the 
classical results must be modified. To this end an arbitrary function f(z, y) must be 
selected which is finite and continuous as are also its first and second derivatives through- 
out D except for the boundary I. It is then possible to write 


F(x, y) = f(x, y) + ua, y). (9) 
From (3), u is a solution of the non-homogeneous equation, 
V'u + ku = —(V*S + kf) = filz, y), (10) 


the second part of which is known, subject to the condition that du/dn = 0 on I. It is 
now possible to develop u in an absolutely and uniformly convergent series in ©; over 


the entire surface: 


u(r, i) = Le®(e, (11) 


where c, are constants in a Fourier series. As u is an analytic function, the series (11) 
can be expressed in a doubly infinite series (a point which is to be verified a posterior). 
Combining (10) and (11), one obtains 


> (kh — ke.®: = fi 


s 
and, with reference to (8), 


e, = Ho bets (12) 
which is useful if the function f, is integrable in the domain D. Thus the solution is seen 
to lead to a unique result. 

Definition of auxiliary function. The relations (9), (11) and (12) resolve completely 
the problem as posed for the non-resonant oscillation in Chap. I of [1] for which Ps Ki, 
with the reservation that f,(z, y) must be defined and the differentiation of (11) term 
by term must be justified. For this problem, the development from (11) and (12), although 
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absolutely convergent throughout the domain, does not have all of the properties of 
regularity. The validity of the formal solution is studied in a subsequent paragraph. 

The definition of the function f, is the next step in the analysis. For a rectangular 
harbor, a complex analytic function X is defined: 


Ne) = 4 [(@ — a) Nn@ — a) — @- 8) n@— 8) = Wile, + ive, »), (13) 


where z = x + iy. The logarithms are evaluated in such a way as to be real if z is real 
and greater than 8. It is evident that \(z) is analytic throughout the domain except for 
the points z = a and z = 8. If the polar notation 


pi = |z—-a|, p2 = |2 — BI, 


6, = arg (z — a), 6, = arg (2 — 8), 


is used, 
dh _ 9¥2_ OW, _ A, wm, tA, _ 
de = ay i] ay = - In ps + pm (0, 62) 
and hence 
f y. 
oy T 


For the portion of r defined by y = 0, the harmonic function y, (z, y) satisfies the following 
conditions: 


dy,/an = 0, O<2r<a,B<2<a) 
dy,/an = A, (a<2< 8B) 


and is regular for the remainder of I. 
Another regular function ¥(z, y) is now sought such that dJ/dn = 0 on the boundary 
y = O and dy/dn = dy,/dn on the other three sides. One can take 


f _ V1 _ y, (14) 
where, from (13), 


y, = a [(z — a) In p, — (x — B) In pz, + y(@ — 4,)). (15) 


. . . . . . - . . 
It is evidently possible to form an infinite number of functions y_which satisfy the 
required conditions; the one chosen herein has certain advantages for numerical calcula- 
tions. For convenience the analytic functions 


dy; dy, 
Noe ee ee ee, Se mae 
p(2) By los ps Os2zSa, 
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are defined. Now, with reference to (14), it can be verified that the function, 
2 2 2 

, = y _ «= i 

W(x, y) = plz) 55 + aw(z z) + r(y) oa 


2 


- { p's + [p'(@) — p')] Z| 8) 


can be chosen. As the order of differentiation is immaterial, g’(0) = r’(0) = 0, q’(b) = 
p’(0), r’(b) = p’(a). 
From (14), (15) and (16), f can be determined, and from (10), 
fi(z,y = V*ytkhy—ky,. (17) 


The quantity V*y + k’y is analytic and regular throughout D and on I. The function 
y, is harmonic and regular except for z = a and z = 8, at which points it is continuous 
but presents singularities of the type 


p, cos 6, In p, — yA, , P2 COS 6, In p2 — yb . 
In the vicinity of z = a, for example, 


he Aina tH, Ma -%o,, 
in which H is an analytic function. 

In the course of the experiments in the Neyrpic Hydraulics Laboratory, vertical 
vortices were observed to form periodically at either end of the entrance (Chap. V of 
[1]). These are the physical counterpart of the mathematical singularities already 
discussed. 

The development of a Fourier series, which is convergent in the domain D, can now 


be written: 


fi(z, y) = > (> Yn,m COS may 7 , (18) 
where 
Yn.m _ I Si Paim de. (19) 
It can be shown that 
c 
,- | < nm , 


in which C is a constant. From (7) and (12) 


— 

nm(n? + m’)’ 

in which C, is a finite number, so that the convergence of (11) is absolute and uniform. 
The development (11) is differentiable term by term, once in x and once in y; the resulting 
series are still absolutely and uniformly convergent in D and I, but the convergence of 
the quantities °u/dz* and 0°u/dy’ is not known to be absolute. 


| Cam |S 
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Resonant motion. From an inspection of (12), the method of the preceding para- 
graph is not valid for the case designated as resonant, for which k’ = k? . It is evident 
that there can be added to the solution for u(z, y) an unknown quantity of the form 
Bo,,,.(x, y) in which B is an arbitrary constant. 

If ds is an element of I’, from (6) and (10) and from Green’s theorem, which is appli- 
cable because of the assumed regularity of wu, 


|| (6, .V u — uV"®,.m) do = i (s... ee u 18.2) ds = /| f.:%,,..d0 = 0. 
Jp r dn dn D 


The existence of the function u to which Green’s theorem is applicable is postulated 
and demonstrated subsequently. 

The necessary condition for the existence of a solution, distinct from B®,,,, in the 
resonant case, is the orthogonality of f, and ®,,,, . This condition is not in general 
fulfilled if A is not zero. From the physical point of view, this would mean that the 
crest of the clapotis must, in general, coincide with the plane of the entrance (so that 
V, = 0) if k = k,,,, . Thus only the resonant solution of the type B@,, ,,(z, y) is possible, 
and the movement is similar to that which would occur for the same harbor shape and 
wave period and with the entrance closed. This result corresponds to that observed in 
the laboratory. Some important though isolated exceptions do, however, occur. 

The partial derivatives of y, are infinite at the singular points but integrable. It can 
be shown that Green’s theorem is applicable to y, and hence to f, throughout the domain 
[3]. From (6), 


fx = -—V'f — Kf 


and, hence, 


|] tiem do == ff nT°f — $T*Fan) do 
7D “D 


B 
/ (s df _ j Han) ds mn | ®, (x, 0) dx 
y 7 


“™ dn “. adn 


2A (a\'’"/ . nxB . Na 
= —|- sin —— — sin —}. 
am \b a a 


Finally, the relationship, 


Il 


. mp . WMNea 
A( sin — — sin — ) = 0 
a a 


is obtained. 
If the dimensions of the harbor do not satisfy one of the conditions 
2N 2N 
B-a_2N B+a_2N+1 it 


™ ’ “ow ’ 


a n n n 


where N is any positive integer, f, and &,,,:can only be orthogonal if A = 0, as already 
stated. If, on the contrary, one of the conditions in.(18) is met, (11) and (12) will give 
the following solution to (10): 


© 


u(x, y) = > c,2; + Be, ; (21) 


‘ 
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where the symbol >~’ indicates that the (n, m) term is excluded from the summation. 
Calculation of the surface displacements. The function f, from (10) was introduced 
only for consideration of the convergence of (11); the singularities are thus isolated, and 


the expressions 
P.@ 
Fy a2, y) = f(a, y) + Qo CnmPnm(2, ¥) 


satisfy for any values of p and q the appropriate boundary conditions for (5). From the 
practical point of view, the values d/,,,/dx and dF,,,/dy are acceptable approximations 
to the values 0F'/dx and dF /dy for small values of p and qg throughout D. The foregoing 
method is thus recommended for the calculation of the velocities from (2). In laboratory 
experiments, however, the point of primary interest is the amplitude of the oscillations 
of the free surface. The relationship from which these are determined is well known: 


_ 14g) 
g Ot |,-0 


From (2), the calculation of z requires only the definition of F(x, y), and not that of 
its derivatives. Development of F in series is hence desirable so as to avoid the difficult 
calculation of the coefficients y,.,, . As F is a regular function, the procedure can be 
developed as follows: 


r 


F= +» an (22) 


n,m 


in which 


i: es I F6, .. de. (23) 


d 


The development of F in series converges absolutely, but, as can be shown, it may 
not be differentiable term by term, and the partial summations in (22) would then not 
satisfy the boundary conditions. From (3) and (23), 


oe [| 6, .V'F de. 
v/D 


Also, Green’s theorem is applicable, so that 


dam = || FU°.mdo + | (Om aF — PSs) ae 
JJp r dn dn 


‘ 2A [a\'?/ . an . wna 
= Kb dhe + (2) sin mn sin — 
an \b a a 


2A a\'*/ . wn . 7 ; 
dn = a (2) (sin Se .. in 2}. 24) 
mk” — ky.) \b a a 


Equations (22), (23) and (24) give the expression sought for ’. The Fourier series 
thus obtained is absolutely and uniformly convergent in D and I, but it is not differenti- 
able term by term. To repeat, the boundary conditions are not satisfied for dF/dn 
calculated by differentiating the development (22). 


or, if k* # kin ; 
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Conclusion. From the equations derived, computations can be made of the local 
amplitudes of both the velocity and the surface displacement. The double-series solution 
obtained for the velocity potential from (9) is absolutely convergent. A more direct 
solution in a double series is provided for the surface displacement. Two classes of 
solutions designated as resonant and non-resonant have been clearly delineated. 

For the non-resonant case, the general solution is given by (9) and (11-16). For the 
resonant motions, the simple solution from the appropriate eigen-function is obtained 
unless one of the conditions in (18) is satisfied. In the latter case (19) is used in place of 
(11). The more direct method for the calculation of the surface displacements is com- 
prised of (22-24). In addition, the basic relationships (1-3) are used in each instance. 

Although the singularities and the lack of continuity for the tangential component 
of velocity at the entrance require further study for the establishment of the analytic 
validity of the solution, the physical validity and usefulness of this type of result were 
demonstrated in [1]. 
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THE PLASTIC INDENTATION OF A LAYER BY A FLAT PUNCH* 


BY 
R. T. SHIELD! 
Fort Halstead, Kent, England 


Summary. Upper and lower bounds for the average pressure in the indentation by 
a flat, smooth punch of the plane surface of a layer of elastic-perfectly plastic material 
resting on a rough, rigid base are obtained by the application of the limit-design theorems. 
The material of the layer is assumed to obey Tresca’s yield criterion of constant maxi- 
mum shearing stress during plastic deformation. The square punch problem is considered 
in detail for layers whose thickness is greater than one-fourteenth of the width of the 
punch. For thinner layers, reasonably close upper and lower bounds for the average 
pressure over the square punch are obtained as functions of the relative thickness of 
the layer. The circular punch is considered briefly, and the bounds obtained determine 
the indentation force with sufficient accuracy for layers which are not too thick compared 
with the width of the punch. 

1. Introduction. The problem under consideration is the plastic indentation by a flat 
smooth punch of the plane surface of a layer of elastic-perfectly plastic material resting 
on a rough rigid base, Fig. 1. A complete solution to the problem under conditions of 


P 


y, // 


/ 


jee © 7 A 


Y ia / 


Fig. 1. Rigid punch on layer of material resting on rough rigid base. 




















plane strain has been obtained by Wang and Lee [1]’. The three dimensional punch 
problem’is one of considerable difficulty at the present stage of elastic-plastic analysis. 
However, as in previous work [2] on the indentation of an infinitely thick layer, the 
limit-design theorems [3, 4, 5] can be used to determine upper and lower bounds for the 
indentation force during the incipient plastic yielding of the layer. 

In the following work the material of the layer is assumed to obey Tresca’s yield 
criterion of constant maximum shearing force k during plastic deformation. This yield 
condition approximates to the behaviour of a frictionless cohesive soil such as clay. 
The square punch is considered first and the bounds obtained for the average pressure 
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on the punch for layers of varying thickness determine the indentation force with 
sufficient accuracy for practical purposes. The lower bound 5k for the punch pressure 
obtained in [2] for a square or rectangular punch on a layer of infinite thickness is 
improved to the plane strain value (2 + r)k. Finally, bounds are obtained for a circular 
punch, and these bounds are in close agreement for layers which are thin compared 
with the width of the punch. The stress field used to obtain lower bounds for the circular 
punch pressure can be adapted to provide lower bounds for any convex area of inden- 
tation. 

2. Limit-design theorems. ‘The limit-design theorems developed recently [3, 4, 5} 
provide a general method for determining the point at which collapse will occur in assem- 
blages of rigid and of elastic-perfectly plastic bodies. Collapse is used here to describe 
conditions under which plastic flow would occur under constant loads if the accompany- 
ing changes in the geometries of the bodies were disregarded. We now outline the limit- 
design theorems in the form required in the following work. 

The boundary conditions are assumed to be of the stress type, i.e. each component, 
T., T,, or T. , of the surface traction is given except where the corresponding velocity 
component, V, , V, , or V, , or the corresponding relative velocity component at the 
interface of an assemblage is prescribed to be zero. A stress field is called statically 
admissible if the field satisfies the stress boundary conditions, satisfies the equations of 
equilibrium throughout the volume of the bodies, and if the yield condition is violated 
nowhere. Where the stress components are continuous functions of position, the equa- 


tions of equilibrium are 


Oc, OT,, OT.; 


+ = 0 (1) 


Ox Oy Oz 
and two similar equations. The equilibrium conditions across a surface of discontinuity 
in the stress field require that the expression 


ON, t+ TryNy t+ 72M, (2) 


and each of the two similar expressions, must have the same value whether it is evaluated 

from the stresses on one or the other side of the discontinuity surface. In expression (2), 

n,,n,,and n, are the components of the unit normal to the surface of stress discontinuity. 
From any given velocity field V, , V, , V, , strain rates can be derived in the usual 

way 

fe] 'g ‘ OV, 


eee ag eee; (3) 
OY 


" Qa 


= —.-- = 
If the strain rates are considered as purely plastic strain rates, the internal rate of dis- 
sipation of energy can be calculated. Also the external rate of dissipation of energy can 
be calculated as the rate at which the applied tractions do work on their points of appli- 
cation. A velocity field is called kinematically admissible if the velocity components 
V., V,, V. satisfy the incompressibility condition 


OV. , OV, i at (4) 


Ox ss Oy Oz 
satisfy any imposed velocity boundary condition, and if the total external rate of dissi- 
pation of energy is not less than the total internal rate of dissipation of energy. 
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The limit-design theorems for the case of all surface tractions increasing in proportion 


have been formulated as follows: 


Theorem 1. Collapse will not occur until the largest values of the surface tractions 
are reached for which it is possible to find a statically admissible stress field. 


Theorem 2. Collapse will occur under the smallest values of the surface tractions 
for which it is possible to find a kinematically admissible velocity field. 

With these two theorems, lower and upper bounds can be obtained for the collapse 
values of the applied surface tractions. 

Since Tresca’s criterion of constant maximum shearing stress is used in the following 
work, statically admissible stress fields must nowhere involve a shearing stress greater 
than the value k. Also, for this material, the internal rate of dissipation of energy due 
to a plastic strain rate is given by [2] as 

D = 2k max | €’ |, (5) 


where max | e’ | denotes the absolute value of the numerically greatest principal com- 
ponent of the plastic strain rate. It now follows from the definition of a kinematically 
admissible velocity field V, , V, , V. , that the velocity field must satisfy the inequality 


[ (7.V.+7,V,+7.V.) dS > [ 2k max |¢' | dV, (6) 
7S “D 

where S is the bounding surface and V the volume of the body or assemblage of bodies, 
and where max | e’ | denotes the absolutely largest principal component of the derived 
plastic strain rate. The strain rate is zero for a rigid body so that rigid bodies in the 
assemblage contribute nothing to the right hand side of the inequality (6). Disconti- 
nuities in the tangential velocity are permissible in kinematically admissible velocity 
fields. For each surface of discontinuity S, in the field, the expression 


[ kavds, (7) 
“SD 

must be added to the right hand side of (6), where AV is the magnitude of the relative 

change in velocity across the surface Sp . 

3. Statement of the problem. We are considering the indentation of the upper plane 
surface of an infinite layer of perfectly plastic material by a flat-ended rigid punch or 
footing. The lower plane surface of the layer is parallel to the upper surface and rests 
on a perfectly rough, rigid base (or equivalently, there is complete attachment between 
the base and material in the layer). The indenting surface of the punch is smooth and 
it is assumed that the free end of the punch is normally and centrally loaded. The force 
P on the punch (Fig. 1) is increased until indentation occurs due to plastic yielding of 
the material in the layer. The punch, the layer and the base constitute an assemblage 
of bodies to which the limit-design theorems are applicable. Statically admissible stress 
fields for this problem must be such that the pressure between the punch and the material 
is normal to the plane surface. Also the resultant normal force must act through the 
centroid of the area of contact, but the distribution of the normal pressure is not restricted 
in any other way, since the material of the punch is rigid. At the interface of the layer 
and the base, the tractions may have any value compatible with statically admissible 
stress fields in the material of the layer, since the base is rigid and is perfectly rough. 
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The rigidity of the punch implies that kinematically admissible velocity fields must 
be such that the area of indentation moves as a plane surface, although there may be 
slip between the material and the punch, which is smooth. The material immediately 
adjacent to the rough base must not move relative to the base. However, a surface of 
velocity discontinuity between the material and the material immediately adjacent to 
the base is permissible. 

4. Square punch. Upper bounds for the indentation pressure. Upper and lower 
bounds for the indentation pressure have been obtained previously [2] for the indentation 
of an infinitely thick layer by a square punch. It was shown that the average pressure 
over the square area of indentation lies between the values 5k and 5.71k. The kinematic- 
ally admissible velocity field used in [2] to obtain the upper bound 5.71k involved non- 
zero velocities in a region of the material of the layer which extended to a depth a sin 34° 
below the surface of the material, where 2a is the width of the punch. It follows that 
5.71k is also an upper bound for the average pressure in the indentation of a layer of 
finite thickness h, provided that the thickness of the layer is such that 


: < cosec 34° = 1.79. 

h 
For layers which are such that a/h is greater than 1.79, the velocity field can be modified 
in an obvious manner so that it lies wholly within the material of the layer and can 
therefore be used to provide upper bounds for the average indentation pressure. The 
relevant formulae are in fact given in [2] so that, in the interests of brevity, only the 
results of carrying out this procedure will be given here. It is found that the values 
5.73k, 5.86k, 6.10k are upper bounds for the average identation pressure when a/h has 
the values 2, 2.5, 3 respectively. These values of the upper bound p, are plotted against 
a/h in Fig. 2, and are joined by a smooth curve, which is marked I in Fig. 2, to the 


line p, = 5.71k(a/h < 1.79). 






































Fic. 2. Upper and lower bounds for the average pressure over the square punch. 


For values of a/h greater than 2.3 approximately, a better upper bound is obtained 
by means of a velocity field which is somewhat similar to the velocity field of [2] men- 
tioned above. The velocity field has two vertical planes of symmetry which pass through 
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the mid-points of opposite sides of the square area of indentation. Hence we need only 
consider that part of the velocity field which accommodates the downward motion of 
one-fourth of the total area of contact. This portion of the velocity field is shown dia- 
grammatically in Fig. 3. The point B in Fig. 3 is a corner of the square area of contact 
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Fic. 3. Velocity field for square punch on layer. 


between the punch and the layer and A, C are the mid-points of the two sides of the 
square which meet at B; O is the centre of the square. The part of the velocity field 
accommodating the downward motion of the area OABC involves flow in the volumes 
BDOEJIHMLKF and BGDPQR, which we shall call volumes 1 and 2 respectively. 
Volume | consists of tetrahedra BDMN, BNML, BNKF, a section of a right circular 
cone BNLK, triangular prisms DMNOHA, NMLAHI, NKFAJE, and a section of a 
right circular cylinder NLKAIJ. Vertical sections of volume 1 by planes parallel to 
the line OA are geometrically similar to the section through OA but vary in size for 
those sections which cut the line BN. The rectangle MLIH is at a distance h below 
the upper surface of the layer of material. Angles AOH, HAI are denoted by a, 6 respec- 
tively, while angles JAH, AEJ are both 45 degrees. The region BDMNOHA moves as 
a rigid body with velocity v cosec a in the direction OH, where v is the initial downward 
velocity of the punch (i.e. the downward velocity of the area OABC). In the region 
BNMLKFAHLJE, the streamlines of the velocity field are parallel to the line HJJE 
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and the velocity along each streamline has the constant value va/h. This value of the 
velocity is obtained from the condition that the relative velocity between the regions 
on either side of the velocity discontinuity surface BMHAN must lie in the plane 
BMHAN. 

In volume 2, the tetrahedron BDPC moves as a rigid body with velocity v cosec 6 
in the direction DP, where 6 is the angle CDP. In the remainder of volume 2, which 
consists of two tetrahedra and a section of a right circular cone, the streamlines of the 
field are parallel to the line DPQRG and the velocity is vc/h, CD and CQ being of lengths 
c and h respectively. 

Energy is dissipated in the discontinuity surfaces between the material at rest and 
the material flowing in volumes 1 and 2, and in the discontinuity surfaces BMHAN, 
BCP. (The rectangle MLIHA is a surface of velocity discontinuity between the material 
in the prism NMLAHTI and the material lying in the lower surface of the layer immedi- 
ately adjacent to the base). Energy is also dissipated in the cylindrical and conical 
regions where the plastic strain rate is non-zero. With expressions (5) and (7), it is not 
difficult to calculate the total rate of internal dissipation of energy. It is found that the rate 
at which work is done by the force on the punch will not be less than the internal rate of 
dissipation of energy if the average pressure on the punch is not less than f(c, 8, y), where 
3 h 


: h 3 ¢ : 1 = 
fc, B, y/k=14+24+2-£-t h(a — 3°) sec! Peal 
fle y) ] aT) h tanB+re > ec B+ 57 se 07 


{ (1+ 5) ( ) 
C 2. ae ne Lad ? 
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2 ,2\1/2 
: +5) (1 + + tan). (8) 


a’ a” 


+ 5 ( 
{ 

In this expression y is the angle PCQ, the angles a, 6 being obtained from the relations 
cot a + tan B = a/h, cot 6 + tany = c/h. It follows from Theorem 2 that the function 
f(c, 8, y) is an upper bound for the collapse value of the average indentation pressure. 
Minimizing the function f(c, 8, y) with respect to c, 8, y provides the best upper bound 
for a particular value of a/h which can be obtained from this velocity field. The following 
table gives the upper bounds p, obtained by approximately minimising f(c, 8, y) for 
certain values of a/h, together with the values of c, 8, y used in each case 


a/h: 2 2.5 3 3.5 4 
p,,/k: De ig 5.81 5.98 6.20 6.47 
c/a: 0.70 0.68 0.62 0.59 0.56 
B: 24° 30° 35° 40° 44° 
v: 19° 23° 26° 27° 29° 


The values of p, are plotted against a/h in Fig. 2 and are joined by a smooth curve 
(marked II in Fig. 2). 

For thinner layers, a/h > 3.4 approximately, better upper bounds are obtained 
with the velocity field shown diagrammatically in Fig. 4. The field is a modification of 
the field of Fig. 3. As before, OABC is one-fourth of the area of indentation and its 
downward motion is accommodated by flow of the material in the regions shown. The 
region which was called volume 2 in describing the velocity field of Fig. 3 is left unchanged, 
while volume 1 now contains four instead of two blocks of material moving as rigid 
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Fig. 4. Modification of velocity field of Fig. 3 to thinner layers. 


bodies under the area OABD. The velocities in the various regions are shown in Fig. 4, 
and are obtained from the velocity conditions at the surface of the punch and at the 
various discontinuity surfaces. It is found that this velocity field will be kinematically 
admissible if the average pressure on the punch is not less than g(c, a, , 8, 8, , 7), where 
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The symbols c, a, , 8, 8; , y have the meanings designated in Fig. 4, while a, , a, a, , 6 
are obtained from the relations a = a, + a, , cot a + tan B = a,/h, cot a, + tan 6, = 
a,/h, cot 6 + tan y = c/h. By Theorem 2, the function g is an upper bound for the 
collapse value of the average indentation pressure, and, for a particular value of a/h, 
the best upper bound is obtained by minimising the function g with respect to the 
variables c, a, , 8, 8; , y. The minimum values of the function g were obtained approxi- 
mately by a process of trial and error for certain values of a/h. The values of the upper 
bound p, obtained and the values of the variables used in each case are shown in the 
following table. It is thought that the exact minimum values of g/k do not differ from 
the values of p,/k shown by more than one unit in the last decimal place. 


a/h: 3.5 4 5 5.5 6 
p,,/k: 6.19 6.34 6.42 6.94 Yes 
c/a: 0.62 0.58 0.54 0.52 0.50 

a,/a: 0.46 0.46 0.48 0.48 0.50 
B: hg 24° 31° 34° 36° 

6. : 25° 28° 33° a5 36° 

y: 28° 30° 34° 34° 36° 


The curve marked III in Fig. 2 joins the values of p, plotted against a/h. 

The velocity field of Fig. 4 was further modified by increasing the number of blocks 
of material under the area OABD from four to six and by increasing the number of blocks 
under the area DBC from two to four. The resulting field provides an upper bound 
which is a function of nine variables. In order to simplify the work, the blocks under 
the areas OABD and DBC were taken to have equal widths in the OA and DC directions 
respectively, and corresponding angles were taken to be equal. The upper bounds 7.18k 
and 7.57k were obtained for the cases a/h = 6 and a/h = 7 respectively. The bounds 
are not the best that can be obtained from this type of velocity field, but it was thought 
that any further lowering of the bounds would not justify the numerical work involved. 
The two bounds are plotted in Fig. 2 and joined by a straight line (marked IV in the 


figure). 
Finally, we describe a velocity field which can be used to give reasonable upper 


bounds for the indentation pressure for large values of a/h. In Fig. 5, LMNO is the 
square area of indentation on the surface of the layer and as before the initial downward 
velocity of this area is taken to be v. The field has four vertical planes of symmetry, 
two passing through the mid-points of opposite sides of the square and two passing 
through opposite corners of the square. Because of this symmetry, it is only necessary 
to consider that part of the velocity field which accommodates the downward motion 
of the triangle CBN. The triangle DEH is the vertical projection of the triangle CBN 
on the lower surface of the layer. BEFGHN is part of a right circular cylinder with 
axis BN, BAFGPN is a right triangular prism, and angle FAB is denoted by a. In region 
CBNHED (and in fact in the whole of the material lying vertically below the square 
LMNO), the velocity components referred to the axes shown in the figure are taken to be 


Pe or an 4 Py eae ee 

ee ee Oe (10) 
These components satisfy the incompressibility condition (4), and satisfy the condition 
that the downward velocity of CBN be v. DEH is a surface of discontinuity between 
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the material flowing in CBNHED and the material at rest immediately adjacent to 
the base. The discontinuity in velocity across this surface is, from (10), given by 


v 2 2\1/2 


The rectangle BEHN is also a surface of velocity discontinuity. The flow in the region 
BAFEHGPN is a plane strain motion with streamlines parallel to the line EFA. The 
velocity has the constant value av/2h along the streamlines in order that the normal 
velocity be continuous across BEHN. The jump in tangential velocity across this surface 


is of amount 


a 2)1/2 
Av = 57 y + 42)". 
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I'ia. 5. Velocity field for wide square punch. 


Energy is dissipated in the discontinuity surfaces DEH, BEHN, HPN, and EFAPGH. 
Energy is also dissipated in the material vertically below CDN and in the cylindrical 
region BFEHGN due to non-zero plastic strain rates. After a little calculation, it is 
found that the velocity field will be kinematically admissible if the average pressure on 
the punch is not less than k times the expression 
; h a 1 so 
ja) = 2 + 2a + cota + Oa (a + cot a) + 0.383 h + (i + 14) 


a h Wy" 2h la ( )"} 
+ 12h log. {2 a + (1 + t n) } + 3 a’ log. is + 1 + 4h? . (11) 
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By Theorem 2, kj(a) is an upper bound for the collapse value of the average pressure 
over the punch. The function j(@) has a minimum value with respect to a when 


(2a + h)' “ 


,* (12 
(4a + h)'” (12) 


sna = 
The upper bounds p, , given by the minimum values of kj(a), are plotted in Fig. 2 
for values of a/h between 5 and 7 (curve marked V in the figure), in order to compare 
the bounds with those obtained above. 

For very large values of a/h, the value of a given by (12) tends to 45 degrees. 
Neglecting the first and higher powers of h/a and taking a = 7/4 in (11), we find that 
for large values of a/h, the upper bound is given by 

- ne a 
Pu = H( 5.07 + 0.383 2), (13) 

Upper bounds for the average collapse pressure under a rectangular punch could 
be obtained by modifying the velocity fields of this section. 

5. Square punch. Lower bounds for the indentation pressure. In this section lower 
bounds are obtained for the value of the average indentation pressure over the square 
area of contact. The statically admissible stress fields employed are extensions into 
three dimensions of two-dimensional stress fields which can be used to give lower bounds 
for the collapse pressure in the corresponding plane strain problem. 

In Fig. 6(a), OH is the centre-line of the punch OA indenting the upper surface of 
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Fic. 6(a). Stress field for the plane strain problem. 


the layer of thickness h under conditions of plane strain. A pressure 5k is produced on 
the surface of the punch by the stress system to the right of the line AF. This stress 
field nowhere violates Tresca’s yield condition and was used previously in [2] to obtain 
the lower bound 5k for the average indentation pressure in the indentation of an infinitely 
thick layer. The value 5k is also a lower bound for a layer of finite thickness. The lines 
AD, AE, and AF in Fig. 6(a) are lines of stress discontinuity and AD, AF are each 
inclined to AE at an angle of 30 degrees. The stresses in the different regions are shown 
in the figure. 

In order to increase the lower bound, the pressure near the centre of the punch is 
increased. This is effected by superimposing a pressure Q in the direction FL on the 
stresses in the region FGKL. FL is then a line of stress discontinuity, and the line KL 
of stress discontinuity separates the stress system in FGKL from the stress system in 
OKL whose principal axes of stress are in the vertical and horizontal directions. We 
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suppose that the regions FGKL and OKL are on the point of yielding in order to increase 
the pressure on OL by as much as possible. Fig. 6(b) shows the Mohr’s circles for regions 
AFL, FGKL, and OKL, which are denoted by a, b, c respectively. The circles a, b for 
regions a, b intersect at the stress point R which gives the stress acting on a line element 
parallel to FL in regions a and b. The pole P, of circle b is the intersection of the line 
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Fic. 6(b). Mohr’s circles for regions (a), (b), (c) in Fig. 6(a). 


joining R to P, , the pole of circle a, with circle b. The stresses acting on line elements 
perpendicular to FL in regions a and b are given by the stress points T and S respec- 
tively, and the geometry of the figure shows that the value of Q is 4k cos 28, where 6 
is the angle GFL in Fig. 6(a). (8 must be less than 45 degrees in order that Q be positive). 
Further, the circles b, c intersect at the stress point S since the pole P, of circle ¢ must 
lie on the g-axis. The line of stress discontinuity KL in Fig. 6(a) is therefore perpendicular 
to the line FL. 

The pressure on OL has the value 5k + 4k cos 28 and the average pressure p over 
OA is found to be given by 


3 = 5 + 4 cos 2341 _ (0.577 + cot a}, (14) 


provided that a/h > 1.577, since 8 < 2/4. The stress field is statically admissible and, 
by Theorem 1, this is a lower bound for the collapse value of the average indentation 
pressure. For a given value of a/h greater than 1.577, the lower bound is found by 
maximising (14) with respect to 8. For values of a/h less than 1.577 the lower bound is 5k. 
For values of a/h greater than 4 approximately, it is better to repeat this process a 
certain number of times, as shown in Fig. 7, so that the pressure distribution on the 
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Fic. 7. Stress field for wide punch (plane strain). 


punch has more than one jump in value. By taking the lines of discontinuity LF, MK, --- 
in Fig. 7 to be inclined at an angle of 224 degrees to the horizontal, the lower bound is 
found to be given by 
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provided that 

a > 
h = 2.991 + 2.828n, 
where n is the number of lines of discontinuity LF, MK, --- . This is a good lower 
bound for large values of a/h. 

Considering now the three-dimensional problem, a statically admissible stress field has 
been given previously [2] for the indentation of an infinitely thick layer by a square punch. 
This stress field shows that 5k is a lower bound for the average collapse pressure and 
the bound also applies to the layer problem considered here. In the material just below 
the area of contact, the stress field involves a vertical compression of amount 5k and an 
all-round horizontal compression of amount 3k. This state of stress extends throughout 
the material in the right pyramid formed by the square area of indentation and four 
triangular faces inclined at 30 degrees to the vertical, the vertex of the pyramid being 
vertically below the centre of the square. In obtaining lower bounds for the layer problem 
we shall assume that the material of the layer which lies in the pyramid is stressed in 
this way. The pressure near the center of the punch will be increased by a method similar 
to that used above in the plane strain problem. 

Since the stress field is symmetrical about the two vertical planes passing through 
the mid-points of opposite sides of the square, we need only consider that part of the 
stress field which supports the pressure on one-quarter of the punch. The stress field 
is also symmetrical about the vertical planes through the diagonals of the square. In 
Fig. 8, OABC is one-quarter of the area of indentation in the upper surface of the layer, 
O being the centre of the square area. For clarity, the points A, B, C, which are shown 
in the plan, are omitted from the perspective in Fig. 8. The triangular faces of the pyramid 
mentioned in the previous paragraph meet the lower surface of the layer in the lines 
FX, RX. The planes LFWT, NRVT are each inclined at an angle 6 to the horizontal, 
and the planes LKPT, NMPT are perpendicular to the planes LFWT, NRVT respec- 
tively. 

As mentioned above, the stress field of [2] involves a vertical compression 5k and an 
all-around horizontal compression 3k in the material in the volume OABCRHFX. For 
the moment we assume that the material vertically above the area HFWVR is stress 
free. An all-around horizontal compression 3k will be added later in this region in order 
to restore equilibrium. In the final stress field the stresses in the vertical plane through 
the line OA will be as shown in Fig. 6(a). [In Fig. 6(a), the line OH intersects the line LK, 
which is not the case in Fig. 8. However, if the width of the punch is such that OH 
intersects LK in Fig. 8, a simple adjustment can be made to the stress field in order that 
the symmetry conditions be satisfied, as they are in Fig. 6(a)]. The minimum principal 
stress in region b of Fig. 6(a) is 3k + 2k cos 26 so that, to satisfy the yield condition, 
the pressure parallel to LT in region LFWTPK of Fig. 8 must have this value at least. 

We assume, therefore, a compression 2k cos 28 parallel to LT in region LFWTPK 
and similarly a compression 2k cos 28 parallel to NT in region NRVTPM. The planes 
TWP, TVP are surfaces of stress discontinuity and we assume that the stress system 
in region TWVP is composed of two equal compressions of amount S in the directions 
TW and TV, so that the plane TVW is stress free. From the equilibrium conditions 
across TWP (and TVP), it is found that S must have the value 


S = 2k cos 28/(1 + sin’ 8) cos’ B, 
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and that a vertical compression of amount 
T = 2k cos 28 sec’ B 


must be added in the material vertically above the planes TWP, TVP. 

Compressions of amount 4k cos 28 are now added in regions LFWTPK and NRVTPM 
parallel to FL and RN respectively. Equilibrium is satisfied across the planes LKPT, 
NMPT by adding a hydrostatic pressure of amount 4k cos 28 in region OLTNMHKP. 
An all-around horizontal compression 3k is now added throughout the material vertically 
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Fic. 8. Stress regions in statically admissible stress field for square punch. 


above the area HF WVR. Vertical compressions of amount 5k in the region vertically 
above the area HFWPVRH and amount 5k — T in the region vertically above the 
area PWUV are added, and it can be verified that the yield condition is not violated 
in these regions. Above the area UV W it is found that the maximum vertical compression 
which can be added without violating yield is given by u(8), where u(8) is the lesser of 
5k and k times the expression 


34 2 _ _2 cos’ 26 sin* 26 - 4 2 cos 28(cos 28 — sin* 8) 
7 (1 + sin’ 8)” cos* B (1 + sin’ 8) cos’ B 








When 8 is less than 25 degrees, u(8) has the value 5k. 
The total force on the area OA BC of the punch due to this stress system is the resultant 
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of a pressure 5k over the area TLABCN less the triangular area vertically above UVW, 
a pressure 5k + 4 cos 26 over the square area OLTN, and a pressure u(8) over the tri- 
angular area vertically above UVW. The average pressure p on the punch is found, 
for a/h > 1.577, to be given by 


> 


cr. ¢ or | h whe | f : Lh 2 afr (2 1 rE 
.@# + 4 cos 26} 1 we (0.577 +. cot 8B)? — oe cot” B{5 — u(B)/k}. (15) 
Since the stress field is statically admissible, this value of p is a lower bound for the 
average indentation pressure, by Theorem 1. For values of a/h less than 2, maximising 
(15) with respect to 8 does not give a lower bound greater than 5k. Lower bounds p, 
for the average collapse pressure for the cases a/h = 2.5, 3, 3.5, 4, 4.5, 5 were obtained 
by maximising (15) with respect to 8. The bounds are plotted against a/h in Fig. 2 
and joined to the line p, = 5k(a/h < 2) by a smooth curve, which is marked (i) in the 
figure. For values of a/h greater than 4.5 approximately, it is better to increase the pres- 
sure near the centre of the punch by repeating this process a suitable number of times, 
as in the plane strain case. Lower bounds for the cases a/h = 5, 6, 7 were calculated by 
repeating the process once and maximising with respect to the angles 6, y of Fig. 7. 
The results are given in Fig. 2, where they are joined by the curve marked (ii). 

When the process of increasing the pressure towards the centre is carried out n times, 


where the n angles 8, y, --- are taken to be 223 degrees, the lower bound is found to be 
la h h? 

Pt = (3.42 +=-+ 1.82 —- — 0.26 ) (16) 
3h a a 

for a/h = 2.991 + 2.828 n. For very large values of a/h, the difference between the 


upper bound (13) and the lower bound (16) is 15 per cent of the lower bound, and the 
true value of the average pressure is given by 0.358ka/h to within +7%. 

We now describe a stress field which provides better lower bounds than those obtained 
above for values of a/h less than 4 approximately. The improvement in the lower bounds 
is less than 2.84% and is therefore not of great interest insofar as numerical results are 
concerned. However, it is thought that the result that, in the indentation of an infinitely 
thick layer, the average indentation pressure over a square punch is not less than the 
average indentation pressure over a long rectangular punch (plane strain), is of sufficient 
interest to justify its inclusion in this paper. The stress field used is an extension into 
three dimensions of a stress field previously used [6] to show the validity of the Prandtl 
stress solution [7] to the plane strain indentation of a layer of infinite thickness. This 
stress field is shown in Fig. 9, in which OD is the centre-line of the punch OA indenting 
the material below the line OF. The region OFED is composed of the Prandtl stress 
field of two regions of constant state, AFE and ADO, and a zone of radial shear AED. 
The Prandtl field is continued to the line of stress discontinuity DSH. The stress field 
below this line is assumed to satisfy the following conditions. (i) ¢, and o, are functions 
of z and z only respectively and r,, is zero, so that the equilibrium equations (1) and 
the symmetry conditions across the centre-line DR are satisfied. (ii) Immediately below 
the line of discontinuity, the material is in a plastic state of stress. These conditions 
and the equilibrium conditions across the line of discontinuity are sufficient to determine 
the stresses in the region DSHRD and to determine the line DSH. It is found that 

—o.(a)/k = 2+ r/2 + 2a — 2 cos (x/4 — a); 
( (17) 


—a,(a)/k = 2/2 + 2a — 2 cos (r/4 — a), 
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where a is the angle GAS, measured clockwise. The angle of inclination y of the line of 
discontinuity to the negative z-axis at the point S is 3r/8 — a/2. Denoting the lengths 
of the lines AD and AS by 1, and r, the line of discontinuity is found to be given by 


r/ro = } cosec” (r/8 + a/2). (18) 
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Fic. 9. Extension of Prandtl stress field for two-dimensional punch. 


As a tends to the value —7/4, the line tends asymptotically to a straight line inclined 
at an angle 7/4 to the negative x-axis. The field is statically admissible and provides 
the lower bound (2 + )k for the average collapse pressure. 

Referring now to Fig. 10(a), O is the centre of the square area of indentation on the 
plane surface of the material, and A, C are the mid-points of the two sides of the square 
which meet at the corner B. The stress field is symmetrical about vertical planes through 
OA and OC and through OB. The pressure on the area of indentation is (2 + )k and 
there is an all-around compression of amount zk in region OABD. The section of Fig. 
10(a) by the vertical plane through OA is given by Fig. 9, and the stresses in this plane 
are the stresses of the field of Fig. 9. This stress system is reproduced, apart from the 
linear scale, in the vertical plane through a line such as LM, parallel to OA. This plane 
is shown in Fig. 10(b). Sections of Fig. 10(a) by planes parallel to the vertical plane 
through OC are similar because of the symmetry about the vertical plane through OB. 
The surface generated by straight lines from the point B to points on the curve DSH is 
a conical surface of stress discontinuity, as is the surface [not shown in Fig. 10(a)] gener- 
ated by straight lines from B to the curve DK. The region vertically below the triangle 
ABD is subjected to stresses a, in the direction OC. These stresses, which are not depend- 
ent on the coordinate y, originate from the discontinuity surface through the curve DK. 
Thus points in the plane ABR, which is inclined at an angle a to the vertical, are acted 
upon by the stress o,(8) given by 


—o,(B)/k = 2 + 2/2 + 28 — 2 cos (x/4 — 8), 
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Fic. 10(b). Vertical section through LM in Fig. 10 (a). 


where is the inclination of CK to the vertical and is given in terms of a by the relation 

cot a = cos B/(2)'” sin’ (r/8 + 8/2), (19) 
KR being parallel to CO. This follows because AR is the orthogonal projection of the 
generator BK (not shown in the figure) on to the plane OAG, which is the plane y = 0. 
As a decreases from 7/4 to zero, the value of 8 given by (19) decreases from 7/4 to 
—7/4. It can be verified that the stress system in the region vertically below the triangle 
ABD does not violate the yield condition. In the same way the region vertically below 
the triangle CBD is in a state of stress which does not violate the yield condition. In 
order to complete the stress field, a compressive stress of amount 2k in the direction AB 
is added to the region x > a, where 2a is the width of the square, and a compression 2k 
in the direction BC is added to the region y < —a. The addition of these two com- 
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pressions ensures that the yield condition is satisfied in the regions x > a, y < —a. 

The field described is statically admissible and it follows from Theorem 1 that 
(2 + m)k, or 5.14k, is a lower bound for the average collapse pressure. For a layer of 
finite thickness, the pressure near the centre of the punch can be increased by the same 
method used in the stress field associated with Fig. 8. The formula corresponding to 


(15) is found to be, for a/h > 2, 


2 2 

: = 5.14 + 4 cos 241 — : (1 + cot »} _ : - cot” y{5 — u(y)/k}, (20) 
where y is the angle corresponding to angle 6 in Fig. 8. Maximising (20) with respect to 
y does not give a lower bound greater than 5.14k for a/h less than 2.5. The bounds 
p, obtained for the cases a/h = 3, 3.5, 4 are plotted in Fig. 2 and joined to the line 
Pi = 5.14k(a/h < 2.5) by a smooth curve marked (iii) in the figure. 
The broken line in Fig. 2 is obtained by averaging the least upper bound and the 
greatest lower bound for the collapse pressure for values of a/h. The upper and lower 
bounds show that this mean value curve determines the collapse pressure to within 
+5% for a/h less than 2, the possible error increasing to +10% for a/h = 7. 

The stress fields used in this section for the square punch can be adjusted in a simple 
manner in order to give lower bounds for the indentation pressure under a rectangular 
punch. 

6. Circular punch. In this section upper and lower bounds are found for the average 
collapse pressure in the indentation of the layer by a flat circular punch. The bounds are 
close for layers which are not too thick in comparison with the width of the punch. 

Fig. 11 shows a simple velocity field which satisfies the boundary conditions on 
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Fig. 11. Velocity field for circular indenter. 


kinematically admissible fields. The field is symmetrical about the z-axis which passes 
through the centre of the punch. In the material vertically below the area of indentation, 
region BCOD in Fig. 11, the radial and vertical components V, , V, of velocity are taken 
to be V, = Vr/2h, V, = —V2/h, where V is the downward velocity of the punch. These 
components represent an incompressible flow. In the remainder of the velocity field, 
region ABD, the streamlines are parallel to the conical surface AD, which is a surface of 
velocity discontinuity. The velocity conditions across the discontinuity surface BD 
together with the incompressibility condition determine the velocity components in 
region ABD to be given by V, = Va’/2rh, V, = Va’ tan a/2rh, where a is the angle 
ADE. A straightforward calculation shows that the velocity field will be kinematically 
admissible if the average pressure on the punch is greater that kq(a), where 


| h ; 7 
eh ee Ee “(1 + 7 tan a) log. (1 + Sent a) -_ i} 
a 3h h a 


+ cosec asec a + tana. (21) 
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Upper bounds p, are obtained by minimising (21) with respect to a for particular values 
of a/h. For a/h = 7, for example, the minimum value of q(a) is near a = 36° where it 
has the value 7.51. If the first and higher powers of h/a are neglected in expression (21) 


for g(a) we obtain 
, la 
ga) = 2+ 3] + cosecaseca + tana. 
€ t 


This expression has a minimum value for a = 35° 16’ so that, for large values of a/h, 


D. = H( 4.83 ‘2 2) (22) 
3h 
is an upper bound for the collapse value of the average identation pressure. 

Attempts were made to extend the stress field of Fig. 6(a) to a circular area of in- 
dentation. It proved to be a matter of some difficulty to provide the appropriate stresses 
normal to the plane of Fig. 6(a) (i.e. the circumferential or hoop stresses) which are 
required by the yield criterion. An alternative two-dimensional stress field, which can 
be extended to the circular punch problem in a simple manner, is shown in Fig. 12(a). 
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Fic. 12(a). Alternative stress field for the plane strain problem. 


The field to the right of the line AF produces a pressure 5k on the portion AG of the 
punch, as in Fig. 6(a). In region FGL the state of stress is a hydrostatic pressure of 
amount 3k, and in region FKL a compressive stress 2k in the direction of FL is added 
to the hydrostatic pressure 3k. The discontinuity line LK is inclined at an angle 6/2 
to the vertical, where 6 is the angle KFL, and a pressure 5k + 2k cos 8 is produced on 
OL. The field is statically admissible and provides lower bounds for the indentation 
pressure. These lower bounds are less than those obtained from the stress field of Fig. 6(a). 

As mentioned in [2], the field to the right of the line AF in Fig. 12(a) can be modified 
so that it applies to any convex area of indentation and to a circular area in particular. 
The modified field lies outside the conical surface through AF, shown in plan in Fig. 
12(b). Vertically above the conical surface through AF, the material is subjected to a 
vertical compression 5k and radial and circumferential compressions 3k. The material 
lying between the cylindrical surface through GF and the conical surface through LK 
is stressed by a hydrostatic pressure 3k. In addition, “legs” of material in this region, 
originating from strip elements of area on the conical surface through LK and inclined 
at an angle 6 to the horizontal, carry compressions of amount 2k. Thus the trapezium 
LKK’'L’ is a strip element of area on the conical surface through LK subtending an 
angle 66 at the axis OH and generating the “leg” LKFF’K’L’ which carries a compression 
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2k in the direction FL. The material inside the conical surface through LK is stressed 
by a vertical compression 5k + 2k cos 6 and equal radial and circumferential compres- 
sions of amount 3k + 2k cos 8. Triangular elements of area, such as L’F’M in the plan, 
on the area of indentation in the annulus bounded by the circles through L and G do 
not lie vertically above the “legs” of material. The pressure over these areas can be in- 
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Fic, 12(b). Plan of stress regions in field for circular punch. 


creased from 3k to 5k by adding a vertical compression 2k in the triangular prisms of 
material below the areas. The field described is statically admissible since it does not 
violate Tresca’s yield criterion and satisfies the equilibrium equations and the boundary 
conditions. It therefore provides a lower bound for the average collapse pressure. For 
layers whose thicknesses are not too large compared with the width of the punch, the 
process of increasing the pressure towards the centre of the punch can be repeated. 


Thus when a/h = 1.577 + 1.414n, where n is an integer, the average pressure over 
the punch is found to be 
je la ~ h? 
pr = KM2.72 +57 + 3.05 —- — 0.44 5], (23) 
3h a a 


where the n angles such as angle 6 in Fig. 12(a) are taken to be 45 degrees. The percentage 
difference between the lower bound (23) and the upper bound (22) tends to zero as 
a/h tends to infinity. 

The stress field used here for a circular punch can be modified so that it applies to 
any area of contact which is convex. It can therefore be used to provide lower bounds 
for any convex area of indentation. 
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ON THE IMPROVEMENT OF ACCURACY IN INTEGRATION* 


BY 
MARK LOTKIN 
Ballistic Research Laboratories, Aberdeen Proving Ground 


1. Introduction. In solving numerically systems of ordinary differential equations 
in a stepwise manner it is frequently advisable, and certainly economically more efficient, 
to vary the step size as the integration proceeds. Further, when applying methods of 
integration that necessitate the computation at the “next” point of a number of approxi- 
mate values for the unknown functions it may be desirable to employ the different 
approximations so as to increase the accuracy obtainable. 

Both of these ideas are not new, and have been used on occasion previously, though 
mostly separately. 2 

The scheme proposed here attempts to combine both techniques in a proper manner. 
The resulting methods are shown to work very effectively on a few critical examples; 
they are also in general sufficiently simple to be easily applicable to modern high speed 
computing machines. 

2. The integration scheme. Let us assume that the integration problem is of the 
form 

y’ = f(x,y), y(%o) = H, (1) 


where y and f(x, y) are column matrices of n components. The associated numerical 
problem is then of the form 


y’ = F(a, Y), 


where Y is a numerical approximation to y, and F(x, Y) differs from f(x, Y) by certain 
rounding errors which are due to the mode of construction of f from its arguments 
x, Y. It will, in general, be convenient to denote mathematically exact quantities by 
small letters, their numerical representatives by corresponding capital letters. Con- 
currently, it should be recognized that the arithmetical operations involving such 
numerical approximants are of the “pseudo-arithmetical” type, of necessity introducing 
certain errors due to unavoidable rounding. 

Let us suppose, further, that the solution of the associated numerical problem has 
proceeded through a sequence of points 2, , X2, +--+ , %-; to a point 2, , so that the last 
step size employed wash = x, — %-,. 

The integration scheme then consists simply in securing two approximations V, , V2 
at %., = x +h, and then utilizing their difference 


éY = V;, od V2 (2) 


for the double purpose of a) improving the accuracy of the V, , and b) deciding whether 
or not a change in step size is warranted. 

The importance of 6Y for accuracy considerations is easily recognized. Let the exact 
value of the solution be denoted by y(x), so that 


Y, = y(u) + ™ ; (3) 


*Received March 23, 1954. 
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m, represents the total error accumulated at x, . The two approximants V, , V, are now 
related to the exact value y(z,,,) as follows: 


Ve = y(n) te, i= | 2 (4) 
To fix the notation let V. denote the more exact determination of y at 2,,; . 


It may be assumed that the exact errors e; are composed linearly of terms due, 
respectively, to truncation, accumulation and rounding: 


e: = > Diatai tt >> Cigtts + > d,,7, . (5) 
a=l1 B=k,k-1 7=1 
Since 6Y = e, — e, it is then clear that a relationship of the form 
ee = abY + >> bite + Dd, chus + > dir, (6) 
a 8 Y 


must exist. The major part of e, , which is easily computed numerically, is 


E,= A8Y + DBT. + YX CeUp. (7) 
e r) 


Let the difference between e, and F, be m.,; . If we now take for the final value of Y at 
2,4, the quantity 
Yin. = V.-8£;, (8) 
then obviously 
Y n+1 = Y(Xp+1) t+ M41 - 

The adjustment of step size may be tied to the magnitude of E, , as follows. Suppose 

it is desired to keep | e, | < ¢, where e > 0 is prescribed. In general, 
| EH, | = me, 

where m is some positive number. It is then reasonable to expect these alternatives: 

(a) if m > 1 shorten h; (b) if m < m) < 1 enlarge h; (c) if m, < m < 1, leaveh 
unchanged. 

Some of the methods for obtaining 6Y may be put conveniently into the following 
two categories: 
I. Application of two distinct methods at full step. The “predictor-corrector” techniques 
[1] are of this type. 
II. Application of a particular method first at full step, then at two successive half 
steps; these may be said to employ “local extrapolation to zero step” [2]. Special cases 
of each of these categories are considered below. 

3. Predictor-corrector method. A good combination is provided by the choice of 


a low order central difference method for predictor, and Heun’s method for corrector. 
The central difference formula is 


Vi = Y,-1 + 2KY{i ; (9) 
its truncation error is t, = (h*/3)y’""(s;), %,-1 < 8: < 2¢s; - Heun’s method, on the other 
hand, employs the formulas: 

Vo=Y¥,+ hVYi, (10) 
V. = Y, + (h/2)(Yi + V3): (11) 
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the truncation error of Eq. (10), which is Euler’s method, is tf) = (h’/2)y’’(8o), that of 


Eq. (11) is tg = — (h?/12)y’""(82), te < 80, 82 < %e+1- It is found that 
é = —h +2pm+umsit-::, (12) 
ep = —t, —2 pth +I +pt+2'pyut:::, (13) 
where p denotes the matrix of elements p;; = df;/dy; , I is the unit matrix of nth order, 
and the omitted terms are due to local rounding only. Since t, ~ — 4 t, , above equations, 
result in 
E, = 6Y — (6T, + PT.) + Ui: ’ (14) 
where 
Une. = (21 + P*)U, — Ui.-1 , (15) 


and the 7; , P are evaluated at 2,4; , V2 


Of interest may be the fact that by subtracting Eq. (13) from Eq. (12) there is 


tp = (1/5)[6Y — (1/2)ptp + (I + p) "uw — w-r]; 


the resulting approximation 


obtained 


tp = 6Y/5 | (16) 


was previously pointed out by Milne [1], p. 28. 

It is clear in general that combinations of formulas having ‘‘comparable” truncation 
errors c,;h"y"" (s;) permit estimates similar to (16). However, combinations not possessing 
this property may still be useful. 

Such a combination is, for example, that composed of formulas (10) and (11). In this 
case one gets 

5 = —t + UI + p)u + oe 
ee = —t+(1+2°'pu + 2'peat+-:-. 
teplacing in the latter equation e, by e. + 6Y there results 
E, = (21 — P)"'[P 6Y¥ — 2T.] + Ui, (17) 
with 
Ure. = (21 — P)'(21 + P/U, . (18) 
Other predictor-corrector combinations may be analyzed in the same vein. 

4. Local extrapolation. For purposes of illustration let us apply the technique of 
local extrapolation to Heun’s method. For the first approximation, then, as by Eqs. 
(10) and (11), 

V, = Y, + (h/2)(Yi + V%). 


Next, let & bisect the interval (2, , x,,,), and let us define 2, = y (&), w = h/2. The 
second approximation V, is then obtained thus: 


W.= Y¥,+ w¥!, 


W, = Y, + (w/2)(Yi + Wd), (19) 
W.=W.itwW, 
V, = W, + (w/2)(Wi + W2). 
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By Eq. (13), where e. now becomes e¢, , 
= —t, —2 pth +I +pt+2 ‘pat: 
Consequently, 
ef = W, — y(&) 
=—7,—-2%'97y+¢(1+r+2 rat, 
where the 7; and 7 differ from the ¢; and p only in that h is replaced by w. Since, therefore, 
T= 2", 7, = 2°*t,, T? i. 
it follows that 
Sel =e, +(71+3p+2'pyua+---. (20) 
On the other hand, a short calculation shows that 
eg = (22 +2 p+2%prei —-(I +2 'p+2"p)u+ -:: 
Therefore, using Eq. (20) while neglecting terms of order p’, there is obtained 
E, = (1/3) + P/3) 6Y + Vis, (21) 
where 
Une = I+ Py. 


It is not difficult to extend this discussion to more accurate integration procedures. 
In doing this let us restrict ourselves to a single differential equation 


y’ = f(x,y), y%) =; 
and utilize the fact that most integration methods are equivalent to forms [3] 


T° = Pe, ¥,,%i,°° TH) (= 1,2,--- WN, 


(22) 
Vo= Yet h¥p+ ++ + UN - DIY, 
V; = x, Go; Y2-; a n(@.V3 + > G,,Y! : i's 

(23) 


+ n'( GeV + > Gy; 1m), 


Here the g,; are constants which in general depend on the range r of points x,_,; to the 
left of x, , and N indicates the degree of the polynomials for which the method is exact. 
Proceeding in the same manner as described above one obtains eventually 


. 14 ge | 2g, + gro) | 
é. = 6Y ———  __ |] 1 + pi te (24) 
’ oN am LT Joo (7”* -l- Joo) (1 + Joo) 
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In those methods where goo = 1, such as Euler’s, Heun’s, Adam’s, Gregory’s, etc., Eq. 


(24) becomes 


bY ont 
ye 1 1 + pg: + 910) oe —F - ***. (25) 


Y—1 


In Heun’s method, for example, 
(Yoo > 91 » 9r0) = (1, 1/2, 1/2), N= 2, 


and consequently 


E, = (1/3)(1 + P/3) 6Y +---, 


which agrees with the principal part of Eq. (21). 
5. Numerical example. The technique described above was tried out on the simple 
case of the differential equation 
y =y’, (26) 


whose solution is y(xz) = (ec — x)~*. For the integration the predictor-corrector combina- 
tion of Euler-Heun was used: 
Vi=Y,4+hYi, Vo = Y, + (h/2)(Yi + Vi), 
E, = (2 — P)“(P 6Y — 27.) + Uns , Ure. = (2 — P)'(2+ PU. , (27) 
Yi4.= V.—E£,. 
For our problem (26), 
P = 2hY., T, = —(h*/2)Y2. 


A. Initial Condition y (0) = 1. 
Starting with y(0) = 1 it was desired to obtain the solution in the intervalO < zr <1 


to an accuracy of three decimal places. 
The results of a number of integrations performed with Heun’s method modified in 


a variety of ways are shown in Tables 1 and 2. 





TABLE 1 
Improvement In Accuracy. Comparisons 
(1) (2) (3) (4) (5) (6) 

x Y (No Corr.) Error Y (One Corr.) Error Y (Two Corr.) Error 

0 1.000000 0.10-6 1.000000 0.10-* 1.000000 0.10-§ 
.04 1.041632 —35 1.041664 —3 1.041667 0 
.08 1.086878 —79 1.086950 —7 1.086957 0 
.12 1.136229 —135 1.136352 —12 1.136363 —1 
16 1.190270 — 206 1.190458 —18 1.190474 —2 
.20 1.249702 —298 1.249972 —28 1.249996 —4 
24 1.315373 —416 1.315749 —40 1.315782 —7 
.28 1.388318 —571 1.388831 —58 1.388875 —14 
.32 1.469815 —773 1.470506 —83 1.470565 —23 
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TABLE 2 
Adjustment In Step Size 








(1) 4 (3) (4) (5) (6) (7) 
x Y 6Y P T> E, E* Error 
0 1.000000 0.10-* 8000.10-5 (— 32.10-§) 0.1076 0.10-¢ 0.1076 
.04 1.041667 — 1632 8333 —38 —35 114 0 
.32 1.470565 — 4406 17635 —150 —125 343 —23 
34 1.515129  —1291 60604 —21 —22 66 —23 
.36 1.562479 —1413 62498 —24 —25 —21 
.60 2.499982 —5528 99994 — 156 —140 469 —18 
61 2.564086 —1582 51281 —22 —23 67 —17 
62 2.631564 — 1708 52631 —24 —25 —15 
.76 4.166679 —6528 83331 —151 —139 454 +12 
765 4.255334 —1827 42553 —20 —22 64 15 
77 1.347844 —1947 43478 —22 —24 18 
. 860 7.148004 —8343 71429 —163 —154 492 147 
. 8625 7.272882  —2298 36364 —22 —24 68 155 
. 8650 7.407570 —2426 37038 —24 —25 163 
. 9000 10.000355 — 5875 50001 —78 —78 244 355 


1. No corrections. For the accuracy required it is clearly sufficient to start with 
h = 0.04. Using Heun’s method without any modifications there are obtained the values 
listed in column (1) of Table 1; the errors, shown in column (2), are seen to exceed the 


tolerance « = 5.10~* already at xz = 0.28. 

2. One corfection. Next the integration was carried out utilizing the expression for 
E, without the term U,,, . The results are exhibited in column (3), the errors in column 
(4). This modification reduces the errors to about one tenth of their previous values. 

3. Two corrections. The inclusion of U,,, in £, produced the results shown in column 
(5); the errors, in column (6), are reduced further considerably. 

The continuation of column (5), together with the appropriate adjustment in step 
size, is shown in Table 2. 

The integration interval was cut whenever the quantity 


E¥ = (2—P)"[|P6Y|+2|/T2/+2+P) (Ui (28) 


exceeded 5.10~*; this occurred at x = 0.36, 0.62, 0.77 and 0.865. At x = 0.9 the step size 
was h = 0.0025, the error 7 = .000355. Use of this step size from the beginning would 
require 360 lines of integration, while the computation shown in Table 2 necessitated 
only about one fifth of that number. 

4. Use of 5Y for step adjustment. Instead of watching £4 for adjusting the step size, 
the quantity 6Y could be used for this purpose, as suggested by Milne [1]. Cutting the 
interval whenever 5Y exceeded 0.0025, which occurred at x = 0.20, 0.48, 0.67, the error 
increased to .000496 already at + = 0.69, so that the use of EH seemed superior, at 
least in this example. 

5. The approximation ZH, = T, . Inspection of columns (4) and (5) shows that 
E, = T, , which is to be expected. If this assumption E, = T, is introduced into Eq. 
(27), there results 
E, = (4 — P)"'P 6bY. (29) 
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An integration based on this formula, with changes of interval whenever | PéY | < 
5.10°*, produced an error of —0.000478 at x = 0.78. A graphical illustration of the 
various techniques is also given in Fig. 1. 
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Further inspection of Table 1 shows that the application of the unmodified Heun 
method leads to an error of —693.10°° at x = 0.60. It is of interest to note that this 
error is quite closely approximated by the error propagation theorem derived in [3]. It 
was shown there that the total accumulated error for Heun’s method is 


k —T, +1; 
Q, = [—1+2°'p, h] Da (30) 


with 
X= A/D = (Goo + pgio)/(1 — pg) = (1 + 2°"p)/(L — 2°"p). 


In our example r; , 72 = 5.10~°, so that the contributions due to these errors may be 


neglected. Consequently, 


= 2 k Kk 2 _ t 
a ~ 2526s) 2G". a 


t=1 





Using in (31) the average value 2~' (p, + p,) for p, there is obtained for the total 
accumulated error at c = 0.60 the value —723.10~°, which differs from the exact error 
—693.10°° by only about 4%. 
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B. Initial Condition y (2) = — 1. The purpose of this computation is to illustrate 
the possibility of continual enlargement of h without the loss of accuracy. The results 
are shown in Table 3. First Heun’s method without any modification was applied, at a 

TABLE 3 
Adjustment In Step Size 


(1) (2) (3) (4) (5) (6) (7) 

z Y (No Corr.) Error Y (Two Corr.) T> E, E* Error 
2.00 — 1.000000 0.10-6 —1.000000 (—32.10~°) 0.10-8 0.10-° 0.10-6 
2.04 — .961568 —30 — .961538 —24 —30 86 0 
2.44 — .694567 — 123 — .694447 —7 —7 25 —3 
2.52 —.658018  —123 — .657893 —48 —57 175 2 
2.60 — .625121 —121 — .625996 —39 —46 1 
3.40 — .416751 —84 — .416664 —8 —8 25 3 
3.56 — .390618 —48 —59 1777 7 
3.72 — .367638 —37 —45 9 
5.16 — . 240379 —7 =—7 25 6 
5.48 — . 223205 —41 —5l 140 9 
5.80 — .208323 —3!l —38 10 
7.72 — . 148801 —8 —9 27 9 
constant step of h = 0.04. The errors are listed in column (2); since now p < 0 the 


errors tend to damp out. 

In the remaining columns (3) to (7) there are again exhibited the computations using 
the technique described in Eqs. (27). The interval was doubled whenever £7? fell below 
25.10°°, which occurred at x = 2.48, 3.48, 5.32, 8.04, --- . The error at x = 7.72 was 
9.10°°, which could hardly have been improved upon even with a constant step size of 
h = 0.04, at a computational labor of 143 lines of integration instead of the 42 needed 


for Table 3. 
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APPLICATION OF THE GALERKIN METHOD 
TO PROBLEMS IN HYDRODYNAMIC STABILITY* 


BY 
RICHARD C. DI PRIMA** 
Massachusetts Institute of Technology 


1. Introduction. Problems in hydrodynamic stability fall essentially into two cate- 
gories: those in which the instability may be traced to the action of viscosity and those 
in which the instability may be attributed to some other factor such as centrifugal force. 
Examples of the former are the nearly parallel, viscous flows, where it is sufficient to 
consider only two dimensional disturbances. In the latter class are the problems of flow 
between rotating, concentric cylinders and over concave walls; here three dimensional 
disturbances must be considered. In these problems there will be no “inner friction 
layer” and hence there should be no rapid oscillation of the eigenfunction. Consequently 
such problems should be amenable to treatment by the Galerkin method. 

The stability of a viscous fluid contained between two concentric rotating cylinders 
was first investigated experimentally and analytically by G. I. Taylor [1]. When the 
cylinders rotate in the same direction investigations have also been carried out by 
Jeffreys [2], Synge [3], Pellew and Southwell [4] and Meksyn [5]. If the cylinders rotate 
in opposite directions, more formidable mathematical problems are presented: however 
the work of G. I. Taylor still applies, though the calculations become much more tedious. 
Also, Meksyn [6, 7] has obtained asymptotic representations of the stability equation 
in this case. 

In Sec. 2, this problem is formulated in a slightly different manner from that adopted 
by Taylor and Meksyn. This formulation brings out clearly the structure of the eigen- 
value as well as its dependence on two other parameters. Explicitly, S, the eigenvalue 
defined in Sec. 2, is a function of Ad, and d,/d. Here 2x/X is the wave length in the axial 
direction, d is the distance between the cylinders and d, is the distance from the inner 
cylinder to the vanishing point of the mean velocity. It is shown, with the help of experi- 
mental data, that for small d,/d, S = S(dd,, 0) and for d,/d greater than roughly .55 
the dependence of S on d,/d cannot be neglected. The Galerkin method is used to solve 
the resultant eigenvalue problem. Computations are carried out for d,/d small, and 
d,/d = 1, obtaining results in agreement with the experimental data given by Taylor 
(see Figs. 1 and 2). The result for small d,/d agrees with that obtained by Meksyn. 
In the transition range from Ad, = constant to Ad = x(obtained when the cylinders 
rotate in the same direction, d,/d > 1) where the dependence on d,/d may not be 
neglected, the method is directly applicable though the calculations will be slightly 
more lengthy. 

The stability of a viscous flow over a concave surface was first investigated by 
Gortler [8] and later by Meksyn [9] who obtained different results. Because of the dis- 
crepancy in the results as well as to see if the Galerkin method might be applicable in 


*Received March 25, 1954. This work was sponsored by the Office of Naval Research under Contract 
N 5 ori-07871. 

**The author wishes to thank Professor C. C. Lin for suggesting this investigation and also for his 
many helpful discussions during the course of the work. Thanks are also due to the Aeroelastic Com- 
puting Laboratory and Miss Ethel Ginsburgh who carried out the computations. 
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Fic. 1. S(Experimental) vs d;/d. 
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a future study of the stability of a compressible flow over curved plates, some rough 
computations were carried out for this problem. The results obtained in Sec. 3 tend to 
confirm those of Gértler. 

2. Rotating cylinders. For a viscous incompressible flow between two coaxial, in- 
finitely long cylinders of radii R, and R, , (R, > R,) and angular velocities 2, and Q, 
the equations of motion admit the stationary solution 








V = Ar+ 2. (2.1) 
r 
The constants A and B, determined by the boundary conditions, are 
a 2 /JQ2 = 2 
A m= a(l - ulee/ Ry) B = ad — (2.2) 
l ae i /R; 1 a R; R; 


where pz = 2,/Q, . 
Superimposing on this basic flow a symmetrical three dimensional disturbance we 
obtain the following linearized differential equation for marginal stability:* 


(DD* — 2) = LN (4 4 ay (2.3) 


Here the wave length of the disturbance in the axial direction is 27/\, v is the kinematic 
viscosity, D = d/dr, D* = D + 1/r, and v is the radial component of the disturbance. 
The boundary conditions are 


v= 0, (DD* — »*)v = 0, D*(DD* — »’)v = 0, (2.4) 
atr = R, andr = R,. 


In accordance with the work of previous investigators it is assumed that the differ- 
ence in the radii of the two cylinders is small compared with their mean, i.e. 


d = (R — R,) K +(R, +- R.). (2.5) 
In this case d is the important scale in the motion, D = 0(1/d), \ = O(1/d) and 1/r 
may be neglected compared with D. Consequently Eqs. (2.3) and (2.4) reduce to the 
following equations given by Jeffreys [2] 


(D? — i’*)*v = AN (4 + 2, (2.6) 


and 
v = 0, (D? — »*jv = 0, D(D*? — »’)v = 0, (2.7) 


atr = R, andr = R,. 

When the cylinders rotate in the same direction the basic angular velocity, A + B/r’, 
may be approximated by a constant. The resultant eigenvalue problem requires the 
minimization of a critical parameter 7'(Q, , u, R, , R.) with respect to Ad[4]. This minimum 
occurs approximately when Ad = z, a result which will be used later. 

For the problem under consideration, in which the cylinders rotate in opposite 
directions, A + B/r’ is approximated by a straight line profile which has the slope of 


*Eliminate u from Eqs. (51) and (54) of [10] after setting H and y = 0. 
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V/r at R, , R, being determined by V(R,) = 0. It is easy to show that 
(1 — »)”" 





Ro = TSR’ (2.8 
to = @ — uRe/R) vated 
Setting r = R, + ¢ we obtain, upon dropping terms of order (¢/R,)’, the relation 
B 2A 
. ee? ee 2 ¢ 
A+ 7 =R. é (2.9) 
Let 
- _ . 8A*d} — 82511 — uR3/R})d? 
£ = oe 1, + ¥ , = r 1 , Ss = a : => — hn nr L —t ae 
1 Zz (d, ¢) a d, PR, Y(1 — R2/RDR, ° (2.10) 
where d, = R, — R, ; then from Eqs. (2.6) and (2.7) 
(D? — a’)*v = Sa*(x — 1), (2.11) 
and 
v= 0, (D? — a’ = 0, D(D* — a’)v = 0, (2.12) 


at x = 0 and z = d/d, > 1, where D = d/dx. 

The stability problem requires, for given values of the physical parameters R, , 
R, , u, the determination of the minimum value of 2, with respect to a. From the form 
of S we see that minimizing Q, for given R, , R, , and » means minimizing S. It is clear 
from Eqs. (2.11) and (2.12) that S is a function of a = dd, , and also d/d, . Hence for 
each new set of physical constants R, , R, , » the minimization with respect to a must 
be carried out anew. Fortunately, however, it is possible for a wide range of values of 
d/d, to neglect the dependence of S on d,/d and then the minimization with respect to 
a may be carried out once and for all. 

The possibility of neglecting the dependence of S on d,/d is motivated by the fact 
that the centrifugal force will tend to cause instability in the region R, < r < R, and 
will tend to stabilize the flow in the region Ry < r < R,.* Now if R, is close to the inner 
cylinder, the exact position of the outer cylinder should have very little effect on the 
eigenvalue problem as long as it is sufficiently far away from R, and hence we can let 
R, — o. That is, for our eigenvalue problem, if d/d, > 1/6 > 1 where 4 is sufficiently 
small we may set d,/d = @ in the boundary conditions (2.12). To check the validity 
of this argument and also to estimate the size of 5, values of S were computed using the 
experimental data of Taylor and plotted against d,/d. The results are exhibited in Fig. 1. 
For 6 somewhere between .5 and .6 we see from Fig. 1 that for d,/d < 6, there is very 
little variation of S with d,/d and the above argument is verified. 

It should be pointed out here that there is possibly an error in the data for R, = 3.80 
({6], p. 491) and probably those points should be slightly higher. 

The analysis, using a linear profile for'V/r, may be summarized in the following 
manner. For 0 < d,/d < 6, S = S(Ad,, 0); for 6 < d,/d < 1, S = S(Ad, , d,/d), and 
for d,/d > 1, T = T(Ad). Then 6 < d,/d < 1 is the transition range from \d, = constant 
to Ad = constant and in that range the dependence on d,/d cannot be neglected. 

First we treat the case in which 0 < d,/d < 6, and the outer boundary conditions 
are to be applied at ~. Then the eigenvalue problem involves only the three boundary 
conditions at x = 0 since we seek solutions which die out exponentially at . This 


*See ({1], p. 327) and [2] for a more detailed discussion. 
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problem is equivalent to that formulated by Meksyn [6]. Formal application of the 
Galerkin method yields the eigenvalue determinant 


| I; _ Sa’ J ;; | = 0, (2.13) 
where 


I;; = [ v(D? — a’)*v; dz, Ji; = [ (x — 1)vyw; dz. (2.14) 

“0 “0 
The functions v,; are taken to be of the form e~™[z* + b(mz’ + 2)] and e"*x" for n > 4, 
where b is determined by the boundary conditions (2.12). 

Because of the change in sign of the term x — 1 in Eq. (2.11) some difficulty in apply- 
ing the Galerkin method is to be expected. In order to approximate the correct mode, 
i.e. the mode corresponding to the smallest positive eigenvalue, using only a few func- 
tions it is necessary to make several preliminary considerations. First from Eq. (2.11) 
we have for x large that v will behave roughly as exp [— (Sa’)'’°x’’*]. Reasonable values 
of S and a from the work of Taylor and Meksyn are 1000 and 2; hence likely values for 
m are 4, 5, or 6. Also taking v; = e "*x" for n > 4 it is easy to show that J;; < 0, and 
cee = 0 as m is 2 n + 1/2; or equivalently as v has its maximum at x = 1 — 1/2m. 
Thus for any single function to give a positive eigenvalue it is necessary that its maxi- 
mum be inside of 1 — 1/2m. This is in accord with the remarks made in the preceding 
paragraphs. 

If all combinations of two functions having maxima within this critical value gave 
satisfactory results, it should be expected that three such functions would give an accu- 
rate value for S, even though extraneous roots might be introduced. For a = 2 the 
following table was constructed 


Function Has maximum at zx = S 
v, =e “[x* + b,(52” + 2)] 0.43 8307 
v, =e [x + b,(42" + 2)] 0.53 3600 
v, = ez" 0.8 7776 
v, =e **z* 0.89 " 
v, =e x 1.0 — 3246 


Two term approximations using various combinations of these v; gave the following 
results: S,,; = 876, S.,; = 968, S., = 1107 where the subscripts refer to the functions 
used. Since the results of these calculations were satisfactory, a three term approxima- 
tion for v was used and the minimum value of S with respect to a was determined 


graphically. The final result was 

Sios = 1075, a= 2.135. 
Meksyn gives S = 1132, a = 2; the mean percentage difference in S is 5 per cent. To 
investigate the effect of using an undesirable function (if only a few are to be used) the 


function v, was replaced by v,; with the result that S,;; = 1200 at a = 2.25. The per- 
centage difference based on S,.; is 11.5 per cent which though not good is acceptable. 
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Both S,.; and the Meksyn results are in agreement with the experimental values of S 
if the data for R, = 3.80, which is probably in error, is not considered. 

For 6 < d,/d < 1 the limiting case when d,/d = 1 was also computed. In this case 
the function (x — 1) does not change sign and no difficulty is encountered in using the 


Galerkin method. Analogously to the work of Pellew and Southwell [4], the function 
P is introduced by 
P = (D° — ay. (2.15) 


The boundary conditions (2.12) are then P = 0, DP = 0,v = Oat xz = Oandz = 1. 
Now, P is expanded in a series of functions P; satisfying P; = DP; = 0 at x = 0 and 
x = 1, and the »p; are obtained by integrating Eq. (2.15) subject to v = 0 at x = 0, 
and x = 1. (Introduction of the function P simplifies considerably the computational 
work. It was not used in the previous case because it was desirable there to have approxi- 
mations for v which did not depend on a.) Using the P; as the weighting functions, we 


obtain the eigenvalue determinant 


| Ii; — Sa°Ji; | = 0, 
where 
Ii; = | {(D? — a’)P,-(D* — a’)P,} dz, 
al al 
J;,2=- | (x — 1){ Dy; Dy; + ayy;} dx + | (Dy;)-y; dx. 
The functions P; are the following: P; = x°(x — 1)* which is symmetric about « = 1/2, 
with maximum there, and P, = 2°(x — 1)*° with minimum at x = 2/5. Carrying out the 


computations we obtain 


S = 3390, a = dd = 3.15, for d,/d = 1. 


This result is within 4.5 per cent of the average value of the experimental data, also 
a = 3.15 is in agreement with the value of Ad obtained by approximating the velocity 
profile by a constant, when the cylinders rotate in the same direction. 

Using this result and the fact that Ad = 2.125 d/d, for 0 < d,/d < 4 it is possible 
to draw a smooth curve for dd vs. d,/d which gives the value of Ad for the interval 
5 < d,/d <1. This curve is shown in Fig. 2. The dotted curve is the extension 
of Ad, = 2.125 past the range of validity of that result, roughly 6 = .55, as observed 
from Fig. 1. The circles indicate the experimental results of Taylor. 

Finally, the method just discussed in which the function P is approximated should 
apply readily for values of d,/d less than one, but still close to one. That is, as long as 
the function x — 1 is negative over the major part of the range. A reasonable lower 
bound for d,/d might be .75. It may be necessary to take an additional term in the approxi- 
mation. 

3. Flow over concave surfaces. The differential equations and boundary conditions 
governing marginal stability of three dimensional disturbances for an incompressible, 


viscous flow over a curved surface are 


(D? — o’)u — (dU/dn)v = 0, 
(3.1) 


o'pUu + (D? — o’)*v = 0, 
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and 

u = 0, v = 0, Dv = 0, (3.2) 
at » = 0, and 7 = © ([8] Eq. (2.7)). In this derivation it was assumed that the com- 
ponent of the basic flow perpendicular to the wall was identically zero. Also the dis- 
turbances were assumed to be independent of x, the arc length along the plate, and of 
wave length \ in the z direction parallel to the plate. In Eqs. (3.1), U is the basic flow; 
5 is the length scale in the y direction; 7 = y/6; 0 = a6 where a = 2r/\; up = 2(U,5/r)? 
(6/r); r is the constant radius of curvature, taken positive for the concave side of the 
surface; and D = d/dp. 

Gortler [8] solves the system of equations (3.1) by using Green’s functions to convert 
them to integral equations, which he then solves approximately. Meksyn [9] eliminates 
v from Eqs. (3.1) and then obtains an asymptotic solution in terms of the large parameters 
o and yw. Actually ¢ = 0(1); this may explain the discrepancy between Meksyn’s results 
and the work of Gértler mentioned in the introduction. 

It is clear from Eqs. (3.1) that the difficulties discussed in the previous section will 
not be present here. However, to use the Galerkin method, it is most convenient to 
retain Eqs. (3.1) as they are and approximate both u and v. Then an n term approxima- 
tion for each of the velocities requires the solution of an eigenvalue determinant of the 
order of 2n, and hence the computational work is more than doubled for each succeeding 
term in the approximations for u and v. 

The basic velocity U is approximated by the straight line profile, U = y forO0 < n < 1, 
U = 1 for » > 1. This was the case considered by Meksyn and as one particular example 
by Gortler. The disturbance velocities u and v are approximated by the simplest functions 
satisfying the appropriate boundary conditions: 


u =e "(Ain + Aan’ + --°), (3.3) 


v=e (Bin + Boy’ + ---). 
The choice of a and b depend upon the number of terms of the series which are to be 
taken. One wants to insure that each individual term dies out sufficiently rapidly as 
n — ©, but also the lower order terms should not die out too rapidly. In this work 
only two terms of each series are used, as primarily the interest was of a qualitative 
nature. The constants were taken to be a = b = 1. 
Applying the Galerkin method in the usual manner the neutral curve given in Fig. 3 





° 01 02 o3 o4 os 06 
aé 


Fie. 3. 








62 RICHARD C. DI PRIMA [Vol. XIII, No. 1 
was obtained. The results are expressed in terms of the momentum thickness 4((9], 
p. 261), R, = U,.6/v. The critical value obtained by Gértler for this case was 
R,(6/r)’” = .83 at a? = .15, by Meksyn R,(6/r)'” = 3.65 at a = .400, and by this 
author R,(6/r)'’? = .39 at a0 = .1. This result* then confirms the magnitude of the 
critical parameter given by Gértler. In performing the present calculations, a second 
critical value of R,(6/r)’”? was obtained for the higher mode. It is interesting to note in 
passing that this value is of the order of 2, i.e., comparable with Meksyn’s value. 
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*After this work was completed, the author received a report by Mr. A. M. O. Smith, On the Growth 
of Taylor-Gértler vortices along highly concave walls, Douglas Aircraft Report No. ES 17110. He makes 
a rather complete investigation of this problem for the Blasius profile, using the Galerkin method with 
different choices of the constants a and b and six term approximations. His result for the-critical param- 
eter is about 2/3 that obtained by Gortler, which again tends to confirm the conclusion given in this 
paper. It should also be mentioned that his equations are more accurate than those used by Gortler. 














THE PLASTIC COLLAPSE OF CYLINDRICAL SHELLS UNDER 
AXIALLY SYMMETRICAL LOADING* 


BY 
E. T. ONAT 
Brown University 


1. Introduction. This paper is concerned with the plastic behavior and the load 
carrying capacities of thin cylindrical shells composed of a plastic-rigid material that 
obeys Tresca’s yield condition and the associated flow rule. The discussion is restricted 
to axially symmetric types of loading and edge support. Under these conditions a yield 
surface is found for the relevant stress resultants. The method employed in the deter- 
mination of the yield surface could be extended to the case of general loading of the 
shells of arbitrary shape. 

Once the yield surface and the associated flow rule are known, the approximations 
from above and below to the load carrying capacity can be obtained by using the theorems 
of limit analysis [1, 2, 3]**. For mathematical simplification it is found convenient to 
approximate the curved yield surface by a polyhedron or other suitable shapes. A simple 
example is treated using this approximate approach. 

2. Notations. Equations of equilibrium. Velocity strains. Let r, y, x be cylindrical 
coordinates, the x axis being horizontal, and let the considered shell be bounded by the 
cylinders r = R + h/2, where R is the radius of the middle surface. An element is cut 
from the shell by two adjacent planes through the zx axis and two adjacent planes per- 
pendicular to the x axis. The stress resultants acting on the sides of the element are 
shown in Fig. 1. 














Fia. 1 


On account of the axial symmetry of loading and support, the membrane shearing 
forces N,, = N,, and the twisting moments M,, = M,, vanish identically and the 
bending moments M, are constant along the circumference. 

As regards the transverse shearing forces, it follows from the assumed rotational 
symmetry that only the forces Q, are different from zero. Furthermore, if the shell 
thickness h is small compared with the radius R, then, the stresses o, and r,, may be 
regarded as small when compared with the stresses ¢, and o, . The state of stress at a 

*Received April 12, 1954. The results presented in this paper were obtained in the course of research 


conducted under Contract DA-19-020-ORD-798. 
**The numbers in the square brackets refer to the Bibliography at the end of the paper. 
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generic point of the shell is essentially two-dimensional and on account of the rotational 
symmetry oc, and o, are the relevant principal stresses. 
These stresses give rise to the principal bending moments 


+h/2 +h/2 


M,=- o,2 dz, M,=- o,2 dz (1) 
—h/2 —h/2 


and to the principal membrane forces 


+h/2 P +h/2 
N, = / o, dz, N, = / a, dz. (2) 
=-h/2 —h/2 
The shearing stress 7,, gives rise to the shear force 
+h/2 
a= - [ Tre Az. (3) 
J —h/2 


If the three-dimensional conditions of equilibrium for local stresses are replaced 
by two-dimensional conditions of equilibrium for the stress resultants, the following 


equations are obtained. 





ins = (0 (or equal to the body force in z direction), 
(4): 
dQ. , Ne _ qdM._o — 
ae +R? “Ge ~%=% 


where 7 is the distributed external load acting in the positive r direction. These equations 
are independent of the mechanical properties of the shell material. 

As regards the incipient plastic flow, following the general theory of thin shells we 
assume that the velocity field is such that the normals of the undeformed middle. surface 
remain normal to the deformed middle surface. Since we are concerned with the axially 
symmetrical deformations of the shell, the component v of the velocity in the circum- 
ferential direction is zero, and we have to consider only the components u and w in the 
z and r directions. 

The expressions for the strain rates at a generic point of the shell situated at a distance 
z from the middle surface are then given by the following equations: 


é, = €& — 2K, &, = €,., Crp = Cre = Czy = O,7 (5) 
where 
Ou w ' 
a= a & =P (6) 
are the principal extension rates of the middle surface and 
A. (7) 


kK => 


ax” 


is the rate of curvature change of the middle surface in the z direction. 
The radial rate of strain follows from the incompressibility of the shell material: 


e, = —(e, + ,). ; (8) 
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3. Plastic behavior of a typical element. Yield condition for stress resultants. The 
perfectly plastic material considered in this paper is supposed to obey Tresca’s yield 
condition. For the case of plane stress which is a good approximation to the present 
problem this yield condition is indicated in Fig. 2. All points on the yield hexagon in 
% Strain-rate vectors 
» f , 

g, 


c B ” 
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this figure represent states of stress for which the maximum shearing stress has the 
constant value of o,/2. 

It is convenient to speak of the set of plastic states represented by interior points 
of the segment AB as the plastic regime AB, etc.; similarly, the plastic state represented 
by the vertex A of the yield hexagon will be the plastic regime A, etc. 

The flow rule associated with this yield condition states that the strain-rate vector 
which has the principal strain-rates e, and e, as components is perpendicular to the sides 
of the yield hexagon for states of stress represented by points other than the vertices 
of the hexagon. At the vertices the strain-rate vector is obtained by the linear combina- 
tion (with positive coefficients) of the strain-rate vectors corresponding to the adjacent 
sides of the hexagon. 

The third principal strain is found from the equation of incompressibility (8). 

Now if the element considered in Fig. 1 is to bend plastically, the yield condition 
and the associated flow rule must be fulfilled in each layer parallel to the middle surface. 

On the other hand Eqs. (5) show that the strain-rate vector in a layer is a function 
of z (indicating the position of the layer) and three parameters—the principal strain 
rates ¢, , «, and the rate of change of curvature of middle surface x. Since the magnitude 
of the strain-rate vector is of no importance in determining the corresponding state of 
stress, there are essentially two parameters; hence e¢,/x, «,/x can be chosen as charac- 
teristic parameters. From the flow rule stated above we conclude that the stresses o, 
and o, are also functions of z and e,/x, €,/k. 

If we then compute the stress-resultants as defined by Eqs. (1) and (2), we see that 
the stress-resultants NV, , N, , M, and M, corresponding to the plastic state of the shell 
element are functions of the two parameters e,/« and ¢«,/x. Therefore if the stress resultants 
N. , N. and M, (which occur in the equilibrium equations) are used as rectangular 
coordinates, the plastic states of the shell element define a surface. This surface will 
be called the yield surface for the stress resultants N, , N, and M, . Leaving to the 
next section the formulation of a boundary value problem with the help of this surface 
and of the equations of equilibrium, we proceed to find the exact shape of the yield 


surface. 
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Consider e,/x and e,/x as rectangular coordinates and assume that x > 0. The vector 
OP with the components ¢,/x and e,/x is the strain-rate vector of the middle surface. 
If on the other hand a vector OQ is constructed as in Fig. 3, it is easily found from the 
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Eqs. (5) that OQ represents the strain-rate vector in the layer situated at a distance z 
from the middle surface of the shell. It is seen from Fig. 3, that 


= = (*— 2) / = age, 
Ce K K 


¢ 


where @ is the angle of inclination of the strain-rate vector. 

Once @ is determined, the corresponding stress state is found from the flow rule; 
for example, in the case considered, 0 < 6 < 2/2, therefore the stress-state corresponding 
to OP (and OQ) is the regime B. Moreover if P is situated in the domain defined by the 
inequalities 

oe £>0 

K “a K 
a similar investigation of strain-rate vectors indicates that every layer of the sheet is 
in the regime B. Similar regions are determined and shown in Fig. 3, the capital letters 
denoting the various plastic regimes existing across the thickness of the shell; for instance 
if P is in the region denoted by C, B, the layers of the shell are either in the plastic 
regime C or B. 

By considering the stress-states corresponding to the domains indicated in Fig. 3 
one can easily obtain the complete yield surface of stress resultants. Previous remarks 
indicate that the region B of Fig. 3 is mapped to one point (B’) of the yield surface, 
where 

N, = N, = agh and M, = 0. 


¢ 
Since B corresponds to a two-parameter family of strain-rate vectors, we must 
expect B’ to be a corner of the yield surface. Similar considerations apply to the domains 
C, D, E, A. 
Inspection shows that interesting information is obtained if the parellelogram shaped 
region in Fig. 3 is considered. In Fig. 4 the part of this region where e,/«x > 0 is indicated. 
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If a generic point P representing the strain in the middle surface is considered the previous 
remarks enable us to determine the stress distribution on the element shown in Fig. 1: 
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where the condition that P should be in the region considered is expressed by the following 
inequalities 
‘i 2 h h 
e<BaqoS4E, ~7<7< 45. (10) 
K K “ 
From (1), (2) and (9) we can obtain the ran stress resultants as functions 


of «,/x and e,/k: 


= 
ooh K ooh 


(11) 
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Equations (11), together with the inequalities (10) define a bounded surface in N, , N, 
and M, space. 

If we put No = ooh, My = ooh’ /4, and eliminate e,/«x and ¢,/x from (11), then we 


obtain: 
ny N, Ne a.) a M, of| ¥: (Ne ad ) | (Ne a = ») 
r(X:, NV. M.) — M,. 1+ -\| N, N. 3 + N, i* 0. (12) 


Moreover the inequalities (10) can be written in terms of stress resultants 


1N 1 N, N N, 
= a ee 2 1 eo | Q. (1: 

2No +2 <No -."°s 3) 

The surface defined by (12) and (13) is bounded by the parabolas B’BC’, C’aD’ and 
D’yB’ (Fig. 6). Three different projections of this surface are also indicated in Figs. 

5a, 5b and 5c. It can be shown that by considering the other half of the parallelogram 
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and also by investigating the case x < 0 one obtains three bounded surfaces which are 
identical with the surface defined by (12) and (13). These surfaces are symmetrically 
situated with respect to M./M, = 0 plane. 

The study of the particular case x = 0 adds four flat pieces bounded by the parabolic 
arcs B’BC’ and C’aD’ and their symmetrical counterparts. Moreover the investigation 
of the case e, = 0 gives two right parabolic cylinders whose bases B’yD and its symmetri- 
cal counterpart are shown in Fig. 5b. These cylinders are tangential to the surfaces 
of the type defined by the Eqs. (12) and (13). With this investigation the determination 
of the yield surface is completed and the resulting convex surface is shown in Figs. 


[5(a, b, c)] and Fig. 6. 
It should be noted that the yield curves obtained by the intersection of the planes 
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N. = 0 and N, = 0 with this yield surface were previously investigated respectively 
by Drucker [4] and Onat and Prager [5]. 

Points within the domain bounded by the yield surface of Fig. 6 represent ‘‘safe’’ 
combinations of axial forces and bending moment, while points on this surface represent 
critical combinations of the stress resultants. 
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Under any such critical combination of N, , N, and M, a certain type of deformation 
becomes possible which is characterized by the ratios between ¢, , ¢, and x. 

If the critical combination is represented by a point (N./No , N,/No ,;>M.z/Mo) on 
the surface defined by Eq. (12), then the corresponding flow mechanism (e,/x, €,/x) 
is obtained from the first two equations of (11): 
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If on the other hand the partial derivatives of F of Eq. (12) are compared with (15), 
one obtains 
: : oF . OF oF 
Noe. : Noep : Mok = =e 2 SO SST 
é. ? Note : Mok = S07 7N,) * aN./No) * X(M./M) Ue 


Equation (16) establishes the following important fact: if a critical combination 
N.,N, , M, is represented by a point H on the boundary of the safe domain, and a 
corresponding plastic deformation by the vector HG with the components Noe, , Noe, 
and M,«, this vector has the direction of the exterior normal to the boundary of the safe 
domain at the point H (Fig. 6). This important result may be interpreted as a special 
form of Hill’s principle of maximum plastic work [6]. 

It is easily shown that this result applies to the other parts of the yield surface. 
However the edges where two analytically different parts of the yield surface meet 
require special attention. If across such an edge the yield surface does not have a con- 
tinuously turning normal, the flow mechanism corresponding to the stress-state repre- 
sented by a point on this edge is not uniquely determined. It can be shown that such 
flow mechanisms are mapped to the domains (C, B), (C, D), (E, F) and (A, F) of e,/« 
and e,/« plane (Fig. 3) and the strain-rate vector with the components Noe, , Noe, and 
Mx is a linear combination (with positive coefficients) of two normals to the adjacent 
surfaces at the considered point. 

In a corner like B’, however, the strain-rate vector is a linear combination (again 
with positive coefficients) of the normals of three surfaces meeting in B’. 

4. The mathematical problem. Examples. Consider a shell of length /, under the 
internally distributed load \-p(x) and the axial force AN per unit length of the shell 
cross-section. 

The basic problem of limit analysis is to determine the value of \ for which the 
plastic collapse of the shell occurs under the various end conditions. 

Physical conditions which must be fulfilled when the plastic collapse occurs require 
the formulation of the following mathematical problem. 

It is required to find the functions N,(x), M,(x), Q.(x) and u(x) and w(x) in such a 
way that 
(1) M, is continuous, with piecewise continuous first derivatives; N, and Q, are piecewise 

continuous; uv and w are continuous and have piecewise continuous first and second 

derivatives. 

(2) N, , M; and Q, satisfy the equations of equilibrium (4), in which N, = AN and p 
must be replaced by Ap. 

(3) The point with coordinates N,/N, , N,/No and M,/M, lies inside or on the yield 
surface of Fig. 6: in the first case dw/dx and du/dx must vanish, in the second case 
the strain-rate vector Noe, , Noe, , Mox defined by the Eqs. (5), (6) and (7) must 
satisfy the normality condition indicated in the previous section. 

(4) The end conditions which will be mentioned later must be satisfied. 


In the problem to be discussed the internal load is uniformly distributed [p(x) = ct = p] 
and the free end of the shell is carrying the axial load N(—1 < N/N, < 0) per unit 
length of the circumference. The other end of the shell is built-in (Fig. 7). First of all 
we assume that N is given, and we shall determine the value of p which causes collapse. 
This is equivalent to the solution of the basic problem. 
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It is obvious from the equilibrium considerations that, when the plastic collapse 
occurs, the active plastic states of the shell elements are situated on the intersection 
of the yield surface with the plane N,/N»o = N/No . Such typical intersections are shown 
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in Fig. 5c. This type of yield curve is essentially a rectangle with the corners rounded 
off by the parabolic ares. The sides of the rectangle are given by the following equations 








N. _ N_ | Pe wy’ | 7 
i 1+ N,’ Sige 1 (x and ¥.* 1, (17) 
where 
-i <=> <6 


As a first approximation we shall neglect the parabolic corners, and consider the 
complete rectangle as the relevant yield curve. The flow rule corresponding to this 
approximate yield curve is such that the vector with the components «x and ¢g is per- 
pendicular to the sides of the rectangle of Fig. 5d. The choice of e, depends on the shape 
of the approximate convex polyhedron. We shall assume the following rule: 


é:€:*«=-—-1:1:0 when N,=N, and -—-M,<M,<M.,, 
(18) 

é.:€,:*=0:-1:0 when N,=-N, and -M,<M,<M.,. 

In the corner where N, = N, , and M, = +M, we specify only the following: 
€:¢:«420:0: 21 (19) 


is a possible flow mechanism. 
Now we can determine the edge conditions more precisely. At the free edge (x = 1) 


M, = Q, = 0, (20) 
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at the built-in edge, w = 0 and either dw/dx = 0 or M, = +M, . In the latter case 
the edge becomes a yield hinge [7]. 
For the problem under consideration let us assume that a hinge develops at the 
built-in edge and 
(M,).<-0 = M,. (21) 
The equilibrium equations (4) and (20), (21) then require that the shell should be in 
the regime QR (Fig. 5d) and therefore 
N, N, 


¢ 


No No 
From (20), (21), (22) and (4) we obtain 
«fh a2) 4% | i (“)'] , 
p=, \lt+ Ni) +? M,| 1 N,) |: (23) 


The corresponding stress-resultant distributions are shown in Fig. 7. 


ar 
No 4 


- Xe 
R 


q 








Fic. 8 


Before these results can be accepted, however, it must be shown that they are asso- 
ciated with the velocity distributions u(x) and w(x) that fulfil the requirements of 


(18) and (19). 
It is obvious that 
ax 
~ QR 


w= az, “u= (24) 
are acceptable velocity distributions and the assumption of the hinge circle at the built-in 
edge is justified. If the second theorem of limit analysis [5] is applied using (24), one 
obtains (23) once again. 
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THE USE OF STABILITY CHARTS IN THE SYNTHESIS 
OF NUMERICAL QUADRATURE FORMULAE* 
BY 
HERBERT M. GURK 


Moore School of Electrical Engineering, University of Pennsylvania 


Abstract. A classification of quadrature formulae is made according to their stability 
chart characteristics. Methods are described for synthesizing formulae with such 
characteristics, emphasis being placed upon those methods used for formulae which 
are applicable in real-time simulation problems. The relations between the stability 
chart characteristics and computed results are also discussed. 

1. Introduction. The notion of the stability chart was first introduced by Gray 
[1] in connection with the problem of digital real-time simulation. He described the 
original use of these charts as a check upon the feasibility of using a particular given 
numerical integration method in solving a system of linear, constant-coefficient differen- 
tial equations. In this paper we are going to describe another use of these charts, the 
development of new quadrature formulae, i.e. formulae to be used in solving differential 
equations. 

2. The concept of the stability chart. The quadrature methods which we consider 
fall in three categories: 

1) The open formula 
M N 
lL, = z. Ajotn-; th Lh Djork-; , 
i=l haha 

where x; means 2(ih) and xf = dx/dt,.;, . Such a formula is denoted as an Oyy 

formula; (or an open formula of type Oy). 
2) The closed or repeated closure formula 


P Q-1 
I, = z. Aj Xn-j + h p bj En-i ’ 


3=1 7=0 


where this is an iterative procedure with the value of xj at any step being obtained 
from 2» of the previous step (at the first step we use an educated guess—some open 
method). This is denoted by rC pg . 

3) The mixed formula OyyCpq is the open method Oyy followed by one application 
of the closed method where all ordinates (x;) in both the open and closed formulae 
use the values computed from the closed and all derivatives (x/) are computed 
using the value computed from the open formula (see Gray [1]). 

In all of the above formulae the a’s and b’s are real. 

Suppose now we are concerned with the problem of integrating an equation of the 
type 

z’ = hz, (1) 

where \ is complex. Assume we have decided upon a fixed quadrature formula of one 

of the types mentioned above. Let h be the interval of integration and assume that we 

are given a set of initial conditions for a function z(t) satisfying an equation of type (1). 

Let z = Ah. Then it is known that there exists a positive integer s, constants c; 


*Received April 13, 1954. The work described here was performed in connection with Contract 
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(i = 1, ---, s) dependent on the initial conditions and complex valued functions r,(z) 
(i = 1, ---, s) such that the computed solution of (1) is 
x(nh) = ms c,[r,(z)]" forn > 0. (2) 
i=] 
For i = 1, ---, 8, let w;(z) be the principal value of log [r,;(z)]. Then the given quadra- 


ture formula yields polynomials P and Q such that 


_ Plexp [—w,()]} 


for: = 1, --- ,8, 
{exp [—w,(z)]} 


» 


where P(0) = 1 and Q(0) = 0 for an open or mixed formula, and P(0) = 1 and Q(0) = 


bo- ¥ 0 for a closed formula [1]. 
*)I[/Q(e-*)] for all complex numbers w. Then the roots of 


Let d(w) = [Pte 
g(w) —-z=0 (3) 
are exactly the values w;(z), 7 = 1, ---, s. We call @ the characteristic function of the 


given quadrature formula. 

Let w,(z) = w,(z), where R[w,(z)] > R[w;(z)] 7 = 1, ---, s. Thus in general w, is a 
many-valued function of z. The stability chart for a given quadrature formula is simply 
a picture of the mapping z — w,(z) in a neighborhood of the origin. 

For the moment let us only consider domains where w, is a single valued function 
of z. In this case we call w,(z), the dominant root of Eq. (3). 

In the problem of real-time simulation our foremost aim is stablity, i.e. 


computed solution 
a se =, 
true solution 


as t—o, 


Because of this and since by (2) x(nh) = c exp [nw,(z)] for large n, we certainly would 
like to have the relation w,(z) = ¢[w,(z)] = z for each z in the domain under considera- 
tion. In particular we want w,(0) = 0, i.e. we certainly want stability for the equation 
x’ = 0. Thus we would like w,(z) = z in a neighborhood of the origin. Looking at this 
from the opposite point of view we would want ¢(w) = w in a neighborhood of the origin 
in the w-plane in which w = w,[¢(w)]. 

The stability chart shows how good an approximation of w,(z) = z is obtained for a 
given quadrature formula. Thus it gives an indication of the stability of the numerical 
solution of a linear differential equation using the given formula. 

3. Synthesis of new formulae. In the synthesis of new formulae our aim has been 
to make w,(z) = z in as large a region around the origin as possible. Thus we want w, 
to be single valued there. We have found, however, that it is much more feasible to 
work with the inverse relations. Thus the main criteria considered in the development 
of new quadrature formulae have been: 

1) w,(0) = 0 

2) ¢(w) = winas large a neighborhood of the origin as possible where w,[¢(w)] = w. 


We first consider 2). 


If ¢(w) = w were to hold exactly in some neighborhood of the origin, then we would 
have 
¢(0) = 0, 

¢'(0) = 1, (4) 


@” (0) = 0 forn > J 
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However, since there are only a finite number of quadrature coefficients which determine 
the coefficients of the polynomial P and Q, we see that the infinite set of conditions (4) 
cannot be fulfilled. 

For an open or closed formula, the Eqs. (4) are linear in the quadrature coefficients. 
If we consider these coefficients as unknown and if there are n of them, then the first n 
linear equations above are independent. Thus they can be solved simultaneously. The 
coefficients so obtained are called “‘classical’’ and the corresponding formula, the “classical 
quadrature formula’’ of the particular type in question. The reason for the name is that 
the formulae so obtained are just those which can be obtained by Taylor series matching. 
For an open or closed formula we call the first k conditions of (4) the first k ‘classical 
conditions”. When a mixed formula is such that the open and closed formulae forming 
it are classical then we also call the mixed formula “classical”. 

Now consider the characteristic function ¢ for a mixed formula. Suppose the formula 
uses n* quadrature coefficients. The coefficients of the polynomials P and Q are functions 
of these n* quadrature coefficients. Let the polynomial coefficients be p; , g; , where 
i= 0,1, ---,r andj = 0,1, ---, s and r and s are the respective degrees of P and Q. 
Let R be the range of all possible (r + s + 2)-tuples (po, pi , -**; De» Go» Qi) ***, Ge) When 
the quadrature coefficients take on all real values. Suppose it is possible to find k < n*, 
such that all points in R can be obtained with k of the quadrature coefficients held 
constant. Then really only N = n* — k quadrature coefficients are necessary to determine 
d. 

Consider the first N of the Eqs. (4) for a mixed formula of a given type with: un- 
determined coefficients. In general all the equations will not be linear; some will be 
quadratic. However, if we can solve them simultaneously and we set the k non-essential 
quadrature coefficients noted above equal to 1, then we get a set of quadrature coefficients 
which we call ‘completely classical’. The corresponding formula is a “completely 
classical’’ mixed formula and indeed the various conditions given by (4) are the “com- 
pletely classical’’ conditions. 

It must be stressed that the completely classical conditions and not the combination 
of open and closed classical conditions must be satisfied by the mixed formula coefficients 
in order that ¢(w) = w. The classical coefficients for a mixed formula will always satisfy 
some completely classical conditions but usually not all (indeed no such case has been 
found). Thus in order to investigate the validity of ¢(w) = w near the origin without 
actually inspecting the stability chart itself one should see if the completely classical 
conditions are satisfied (or almost satisfied). 

Another important point is that a completely classical mixed formula may not exist 
for every method OyyCrs . And indeed when one does exist there may appear to be 
more than one such formula (since our simultaneous equations are quadratic). It is 
necessary that for each formula found, classical, completely classical or otherwise, the 
original condition be satisfied, i.e. the condition that w,(0) = 0. Although necessary and 
sufficient conditions that this be satisfied have not been found, several necessary con- 
ditions have been found and may be used as a quick check on the possibility of a stability 
chart corresponding to a useful quadrature formula. We might list a few of the simpler 
ones: 


If 


N 
7*P(1) = R(r) = r* — DY ar 


i=l 
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then the following must hold: 


I. ay | <1 unless R(r) = (r — 1)", 
adR(r : 
PE RW) =N- Je (N — j)a; > 0, 
Gr jaa —_ 
IR) y ] >0 WN odd 
III. — = N(-1)""' — SO (N — pa(-1)*" 
a ~— <0 WN even 


Other similar conditions can also be found. 

Obviously any one of the classical open or closed or completely classical mixed 
formulae corresponds to the stability chart which is most nearly the identity map at 
the origin among the charts for formulae of the same type. However, as mentioned 
above, the condition w,(0) = 0 is not always satisfied for such formulae, and even when 
it is, the neighborhood of the origin where w,(z) = z may be quite small. Thus we would 
want to develop other formulae which would not suffer from these deficiencies. 

As noted previously the only way we can be sure that a formula does not have the 
first deficiency is to check directly, i.e. we must know all the coefficients. For the second 
problem, however, several methods have been used with a reasonable amount of success. 
These are the pseudo-classical or pseudo-completely classical methods (from now on we 
only use the term classical even though the same procedures are also used in completely 
classical cases) and the shifting technique. 

The pseudo-classical method consists of imposing some but not all of the classical 
conditions and then imposing enough other conditions to get the necessary number of 
equations in order to obtain all the quadrature coefficients. The classical conditions 
are used to keep ¢ under control at the origin, the more such conditions taken, of course, 
the better ¢(w) = w close to the origin. The other conditions are used to help control 
¢ at points other than the origin. It is these that we must select carefully in order to get 
a “good”’ quadrature formula. 

There are many tools that may be used in making the choice of the non-classical 
conditions in these pseudo-classical methods. Recall that the branch contour of a 
stability chart is the locus of points z,, in the z-plane such that there exist complex 
numbers w; , Ww. , W; * w, such that ¢(w,) = ¢(w.) = z, and Rw, = Rw, = Rw,(z,.), 
i.e. points where w, is not single valued [1]. In order to make the region where w,(z) = z 
have a reasonable size, it would be advantageous if it could be shown that the branch 
contour is not too near the origin. This suggests the use of the theory of schlicht func- 
tions. Indeed it was found that the wealth of material about regions in which a function 
is schlicht (e.g. [2]) was quite useful in putting bounds on the coefficients remaining 
after the classical conditions had been used to eliminate as many as possible. Some 
bounds have also been found directly by algebraic methods, by determination of the 
end points of the branch contour. 

The criteria used in these pseudo-classical methods may also depend on the particular 
problem for which the quadrature formulae are being developed. For example, in the 
original problem for which this method was used (real-time simulation of a certain 
phenomenon), we were interested only in values of z = x + iy where x < 0. For some 
types of quadrature formulae it was found that we could not make ¢(w) = w throughout 
our whole region of interest. Therefore we tried to increase the size of the region where 
w= u+iv,u < Ocorresponds toaz = x + ty, x < 0. Thus we were assured of damping 
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when there should be some (even if the true and computed rates of damping were quite 
different). Bounds were found on certain of the non-classical parameters (i.e. those 
remaining after classical conditions are used). The bounds themselves represented values 
such that the intersection of the u = 0 and x = 0 lines was a maximum (in y) or such 
that (77) (the intersection of the u = 0 and y = 0 lines) was a minimum (in x)—this 
was a desirable property in open or mixed type formulae. These were of particular 
interest in our problem. In other problems, there might be other aims. In general, 
bounds might be found for which values a particular type of formula might be as nearly 
classical as possible, when if cannot be strictly so and still satisfy the condition of 
w4(0) = 0. 

It must be emphasized that such bounds must be found separately for each type of 
quadrature formula. Indeed, each pseudo-classical method must be tried separately for 
different types of formulae. So far we have not been able to use any one method for a 
general type of formula. 

The only-method of synthesis that is not pseudo-classical and has met with success 
is the shifting technique. It has been observed that in many cases stability charts for 
various formulae show a good approximation to the identity map for some portion of 
the z-plane to the right of the imaginary axis (e.g. classical O,.C,. — Fig. 1). In some 
cases this ‘“‘good”’ first quadrant region shows a better approximation than the “‘good”’ 
second quadrant region in which we were interested in the real time simulation problem. 
It may also be the same for regions of interest for other problems. Thus we are led to 
the following type of translation: 

Let wy) be real such that ¢(w) = w is very good in a neighborhood of w, . 


Let w* = w— WH, O(w*) = o(w) — d(w). 


We want @ to be the characteristic function of a usable quadrature formula. This 
will be so provided that 
1) a quadrature formula having @ as characteristic function can be found; 
2) the dominant root (i.e. w with largest real part) of the equation 6(w) = 0 is w = 0. 
As noted before, condition 2) must be checked in any case. Thus condition 1) is all 
we need worry about. But 1) can always be satisfied since 





P P(e") 
O(w*) = d(w) — o(w) = —— — (wp) 
Ss ake 
_ Ple*” - e**) — p(w) Qe" - e**) _ P*(e“*’) 
Qle-** - e**) Q*(c-”’) ’ 


where P*, Q* are polynomials with real coefficients (this was the reason for w,) being 
real) of degree V, M respectively, and 


N 
N 


M = deg Q, if deg P < deg Q, 
deg P, M = deg Q, if deg P > deg Q. 


II 


Also P*(0) can be chosen equal to 1 unless 1 — ¢(w,.) Q(0) = 0; which can be avoided 
by a slight change in w, . 


'bo = 0 if ¢(w) corresponds to an open or mixed formula 


Q*(0) =), Uo; 
\bo = 0 if ¢(w) corresponds to a closed formula. 
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Fig. 1. Stability chart for classical 01.C 12 


Thus if 
P*(r) = 1+ > a;r", Q*(r) = = 6,7’, 


i=l 7=0 


then 6 is the characteristic function of 
N M 
tr, = >, (—a)a,-; th > d,xh_; , 
i=1 7=6 


which is really a quadrature formula since the a; and b; are real. Such a formula will 
be useful provided N and M are not too large. 

4. Remarks. In conclusion we wish to reiterate that the new formulae were developed 
for a special purpose—real time simulation—where stability is the essential criterion. 
However, for any quadrature problem, stability is still very important even if it is not 
the first criterion to be considered. Thus it seems logical that the methods of synthesis 
outlined above should be combined with other methods in the development of usable 
formulae for any such problem. 
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ON THE EVALUATION OF INDEFINITE INTEGRALS INVOLVING THE 
SPECIAL FUNCTIONS: DEVELOPMENT OF METHOD* 


BY 


L. C. MAXIMON anp G. W. MORGAN 


Brown University 


Summary. ‘The solution of a general second order linear partial differential equation 
in two variables by Laplace transforms is utilized to develop a method for the evaluation 
of a large class of indefinite integrals involving the special functions. 

1. Introduction. Indefinite integrals involving the special functions are of im- 
portance in analysis and in many physical problems. Although extensive work has 
been done on integrals containing the Bessel functions the literature contains relatively 
few results for any of the other special functions. For references see the paper immediately 
succeeding this one in this issue.’ In the following we develop a method of evaluating 
integrals which is equally applicable to many of the special functions, including the 
Bessel functions as a particular case. 

2. Motivation. Let us consider the solution of a linear homogeneous partial differen- 
tial equation in one space variable and time which is typical of those encountered in 


mathematical physics: 


dy _ @ dy 
dt aa | a au. (1) 


We may solve this equation by use of the Laplace transform defined by 
Y(z,9) = [ e-"y(x, t) dt (2) 
0 
in terms of which Eq. (1) may be written 
a(x) Y"’ + a(x) Y’ — sY = —y(z, 0), (3) 
where primes denote derivatives with respect to x. The ordinary inhomogeneous differ- 
ential equation (3) may be solved by variation of parameters, giving 
. Y,y(z, 0 


Y,y(z, 0) —— Yiy(x 
'y a(xW(a) 2 J a(x) W(2) 


in which Y, and Y, are any two linearly independent solutions of the homogeneous 


Y(z, 8) = } 


equation 
a(x) Y’"’ + a(x) Y’ —sY =0 (5) 
and 
W(x) = Y,Y;-— YiY; . 6 
is the Wronskian of Y, and Y, . Thus, assuming that the boundary and initial conditions 
pertaining to Eq. (3) are known, the limits on the integrals in Eq. (4) may be determined 


*Received April 20, 1954. 
1],. C. Maximon, On the evaluation of indefinite integrals involving the special functions: application of 
method, quoted as [B] in the following. 
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and the transform Y(z, s) of the solution is then given in terms of indefinite integrals 


of known functions. 

Rather than proceeding to the inversion of Y(x, s), however, we shall consider the 
integrals in Eq. (4) as our unknowns and hence would like to have Y(z, s) a known 
solution. One case in which the solution is quite simple is provided by assuming the 
boundary and initial conditions to be such that the solution y(z, ¢) is actually time 


independent, so that 


y (x, t) = y(x, 0) = Yo; 


Y(z,s) = “* 


and yp satisfies 
[a(x)y’]’ = 0. (8) 


The particular boundary conditions that must be chosen need not concern us here, 
since, as will be seen later, they do not enter the fiaal formulas. Writing 


I; =8 | zat dx ¢ = I, 2, (9) 
we now solve for J; from Eqs. (4) and (7): 
YiJ.— Yel; = yo. (10a) 
Differentiation gives 
YI, — Yo, = yo . (10b) 
Hence 
I, = ¥Y: -. YoYi 519 (11) 


and from Eqs. (9) and (11) we have 
s | Yiyodx = a(2)[y¥! — ys¥i] i=1,2 (12) 


afver substituting 


= . a’(zx) const. 
W = const. X exp| — | [> |  —-— 
const. X exp | a(z) ax an) (13) 
from Eq. (3). By differentiating Eq. (12) it may be noted that this equation is valid 
for any solutions Y; and y, of Eqs. (5) and (8) respectively. 
Let us consider a particular example. With 
a(z) = 1-2’, s = —v(v + 1) (14) 
we may set 


the Legendre function, and from Eq. (8) 


ss} (16) 


o=A 3 ae = 
y + B log | Ts | 
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where A and B are arbitrary constants. Substituting Eqs. (14), (15) and (16) in Eq. 
(12) we obtain two indefinite integrals involving the Legendre function. For A = 1, 
B = 0, we have 


eee il P(x) de = (1 — 2)P(a), (17) 


and for A = 0, B = | 








. | at 1 
-y(vy + 1) | P,(x) log t =i dx = (1 - 2) Pra log LU = =| = $PAa) I. (18) 
‘ l—2Zz j}l—2 1 
3. Derivation of generalized formula. We now extend Eq. (12) by considering a more 
general equation than Eq. (1), and hence let y(z, ¢) be a solution of 


; : 9 | ay _ 
} Ee + b;(x) | ro fo 0, (19) 


1=0 





in which a,(x) and b,(x) are arbitrary. Introducing higher derivatives with respect to ¢ 
does not seem to generalize further the final results. Additional results may be obtained, 
however, by considering an equation containing higher derivatives with respect to 
x and this is done in [B]. Defining Y(z, s) as in Eq. (2) and taking the Laplace transform 
of Eq. (19) with respect to ¢ we obtain 


> fa,(x) + sb,( (ay) * -> b,(: 7) Se K(z), (20) 
which we solve by variation of parameters, giving 
| Y:K@)_ _ y f__YK@) ' 
Me, = Ms / (a> + sb.) W Ms (a, + + sb.)W dz. (21) 


Here Y, and Y, are any two linearly independent solutions of the homogeneous equation 


y 


> fax) + sb,(x)] aY = (0 (22) 


i=0 


and W is the Wronskian of Y, and Y, , defined in Eq. (6). We again assume that the 
boundary and initial conditions are such that y(z, #) is independent of time, so that 
Eq. (7) is valid and yp is a solution of 


S o{o £8 = 0. (23) 


—_ dz* 


Substituting Eq. (7) in Eq. (21) we have 


Tas Y,/; oa Y,/, ’ (24) 
where 
tae (eae 61,2, 25a) 
CW 


c; = a(x) + sb,(x) ¢ = 0, 1, 2. (25b) 
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As in Eqs. (10), (11) and (12) we differentiate Eq. (24), solve for J; , and equate the 
result to the right hand side of Eq. (25a). In this equation we may then substitute 


W = const. X exp E [() az | (26) 


and 


sK(z) = > (c; — a;) os (27) 


i=0 


from Eqs. (20) and (25b), obtaining finally 
y c, \ l : d'Yo , . | (::) | 
. 1} ly: 2 ‘ i. , ; * , 2X L) dx 
| Y, exp i (<: d. | ra d (c; — a;) ae! dx = (y¥sY; — yoY1) exp | = la 
(28) 


and a similar and redundant equation with Y, in place of Y, . In Eq. (28) Y, is a solution 
of 


2 4, Se = @. (29) 
‘<0 dx; 
Also, yo is a solution of Eq. (23), in view of which Eq. (28) may also be written in the 
form 
/ Y, exp il (2) d | is > d'Yo dz = (YoY, — YoY:) ex il (::) | (30) 
1 exp | / oe Ot Gg SE = (YO, — Yok: exp oa) O° : 


Moreover, in Eq. (30) y may be completely arbitrary; Eq. (23) need no longer be 
assumed. Frequently, however, it will be found more convenient to use Eq. (28), in 
which case Eq. (23) must be satisfied. By differentiating it may be noted that Eq. (28) 
is valid for any solutions Y, and yp of Eqs. (29) and (23) respectively, and Eq. (30) is 
valid for any solution Y, of Eq. (29) and arbitrary y, . Equations (28) and (30) are the 
desired generalizations of Eq. (12). Their use will be illustrated later by some examples. 
We are indebted to Dr. Paul Germain for pointing out to us that Eq. (30) may 
also be obtained by considering the Wronskian 
W* = Yiyis — Yin, 
where Y, is a solution of Eq. (29) and Y, is arbitrary. Thus, if we let 
f(x) -— CoY a’ + C1Yo + CoYo » 
then 
c,W*’ +¢,W* =fY,, 


and solving for W* we have 


W* exp | (:) az | = | man exp | (<1) ax | dz, 
2 2 2 


which is identical with Eq. (30). Although this derivation is far simpler than that given 
in arriving at Eqs. (28) and (30), the lengthier discussion is presented for a number of 


os 
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reasons. First, it is in Eq. (28) that one sees clearly that to obtain fairly simple integrands 
(and hence useful integrals) one must choose Y, and yp so that all but one of the factors 
c; — a, is zero, and this equation is only put in evidence by the first method of derivation. 
Moreover, although relations involving the Wronskian have often been used in the past 
to evaluate indefinite integrals, the exhibition of Eq. (30) solely as an expression for the 
Wronskian W* does not point up the possibility of immediately adapting this equation 
to the evaluation of a very general class of useful integrals. Second, the idea of starting 
with a relatively difficult time dependent problem leading to a solution in terms of inte- 
grals which must be calculated, and then looking for a particularly simple solution 
(e.g. a time independent solution) which may be arrived at by other means, thus enabling 
us to reverse the usual role of known and unknown quantities in the time dependent 
solution, may prove fruitful in connection with other problems. 

While extensive work has been done on indefinite integrals involving the Bessel 
functions, many of which will be obtained by the methods developed here in [B], we 
give two examples of other integrals to illustrate the simplicity of the present procedure.” 


Thus if we choose Y, = W,., , the Whittaker function which is a solution of 
yr 4(-145 +4 =r = 0, 
and if we let yo = W,., , then from Eq. (28) we have immediately 
ot an a il x°W,..W.., dz = W!.,W.., — W.W!, - (31) 
in the other hand, if we let yo = W,., , then again from Eq. (28) we have 
(x — dA) [ « 'W...W., dz = Wi_.W... — Wi... - (32) 


(31) and (32) were arrived at by choosing Y, and y so that 

. (ec; — a,) d'y,/dx’ would reduce to a single term. Evidently many integrals may be 
ev aluate d in this way, but the class of integrals which may be obtained by this procedure 
may be greatly extended by allowing more than one of the terms c; — a; to be non-zero. 
We select functions Y, and y) containing one or more parameters so that they have the 
same functional form but differ in the value of their parameters. Using the recurrence 
formulas satisfied by these functions one may obtain recurrence formulas for integrals; 
the first few integrals may be evaluated in the manner illustrated by Eqs. (31) and (32) 
and the recurrence relations then give all others. This procedure is developed in [B], 
in which a number of integrals involving the Legendre functions, the confluent hyper- 
ge ometric func tion and the Bessel functions are evaluated. 


*See Re fe. [7] and [8] at the end of [B]. 
Vote added in proof: Work of a similar nature by G. M. Muller has just appeared in Proc. Amer. 


Math. Soc., Oct. 1954, p. 716. 
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ON THE EVALUATION OF INDEFINITE INTEGRALS INVOLVING THE 
SPECIAL FUNCTIONS: APPLICATION OF METHOD* 
BY 
L. C. MAXIMON 


Brown University 


I. Introduction. As indicated at the end of the paper immediately preceding this one 
in this issue’, the procedure developed there for the evaluation of indefinite integrals 
will be extended here by permitting a more flexible choice of the functions Y, and yo 
and by making use of the recurrence relations satisfied by these functions to derive 
recurrence relations for integrals, from which these integrals may then be evaluated. 

II. Associated Legendre functions. Let us choose Y, = P%,(x), a solution of Le- 


gendre’s equation: . 
(1 —2)Y" — 2rY’ + Ee +1) - e ly = 0 (1) 
mite 


and let yo = «*P*,(x). Then from Eq. (30) of [A] we have, making use of the differential 


equation satisfied by P*?,(x), 
[ P1200 — a )x'(P*.)’ + MA — Da ?PS, + [vn + 1 


— p(y. + 1) — AA + Da’ P*,} de = (1 — 2’) [P* (2*P*.)’ — a P* (P% )’]. 
Using the recurrence relation 


(—a2)(P)’ = -wPi+vt+wP-s; (3) 


Eq. (2) may be written in the form of a recurrence relation between integrals: 


ih + ~&m® +A +h Ss »] | 2*P* P*. dx 


+ 2rX(v. + ») | a Pt Pk, dx + MA — | z**p* P*, dz 


(4) 
= 2'[(v, — v2 — A)aPY.P?, + (vo + w)P?Pr-1 — it w)PrPr, 
+ r\2"P? P*,]. 
In Eq. (4) we may let A = 0, 1, 2, --- successively and thus obtain 
| STF 2, 48 for n=0,1,2, --- (5) 


provided that none of the denominators in the resulting expressions vanishes. This 


can be shown to be the case for arbitrary wu if 


Vi — wm ¥N, n — 2, n—4,--- —(n— 2), —n 
(6) 
vi tw #~n— 1, n—-3,°--- — (n — 1), -—(n+ 1). 
*Received April 20, 1954. 


1,, C. Maximon and G. W. Morgan, On the evaluation of indefinite integrals involving the special 


functions: development of method, quoted as [A] in the following. 
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Seven additional recurrence formulas for the integrals (5), each also requiring condition 
(6), are obtained directly from Eq. (4), from the possible combinations of (a) inter- 
changing the subscripts v, and v, in the Legendre functions in Eq. (4), and (b) substituting 


—v; — 1 for v; and noting that P*, = P*,,-, , 7 = 1, 2. It may be noted that if any 
one of the numbers vy, , —», — 1, » , —v2 — 1, is equal to — yz, then in two of these 
formulas the second integral in Eq. (4) will not appear and a particularly simple recur- 
rence formula follows. In this case, if —u = v, or —v, — 1 then none of the denominators 
W il] vanish, provided 

utn, p—vw—1#n, n—2,°---, 
and if —u = », or —v, — 1, none of the denominators vanishes provided 

h+n, p—m—l xn, R= 2, 27% 5 


the above series terminating in each case with 0 or 1 depending on whether n is even 


or odd. 
. . . ~ . . h 
A more symmetric recurrence relation for the integrals (5), involving only z° and 
x" in the integrand, may be obtained by interchanging », and vy, in Eq. (4) and sub- 
tracting the result from Eq. (4), giving 


[v.i(vy, +A+ 1) - vv. +rA+ 1)] | 2’P* P*. dx 


‘ 


+ (vo + p) [x 'P? P* _, dx — Mn, + w) | zg" Ps Pe ., dz (7) 


= x [x(y, — v2)P%,Pr, + (v2 + w)PrPr.-1 — (1 + u)P2,P?,-1)- 


A similar equation follows upon substituting —»; — 1 for »; and noting P), = 
p*,,-, (¢ = 1, 2), in each term in Eq. (7). In this latter equation as well as in Eq. (7), 


however, the condition that the denominators do not vanish is, for arbitrary yu, 


Vi —% #N, n-—l1, n—2,-*+—mn 


(8) 
vy tw #~n—1, n—3,--- —-(n—- 1), -m+ 1) 


rather than (6). 
III. Confluent hypergeometric function. As a second example of the type of integral 
which may be evaluated by this procedure we may choose Y, = F(a, ; c; x), a solution of 


xrY¥Y’+(—2)Y’ -—a,Y=0 (9) 
and let yo = «F(a, ; c; x). Then from Eq. (30) of [A] we have, using the differential 
equation satisfied by F(a, ; c; x), 

[ e-*a°~'F(a, 3c; 2)[2dAa° F(a, 30; 2) + MA + — 1)2"' Fae 5 €; 2) 

, + (a, — a, — A)a*F(a, 3c; x)] dx (10) 
= e*2°{F(a, ;c; 2)[2*F(a. 3c; 2)]’ — 2’F(az 5c; 2) F(a, 3c; 2}. 

Using the recurrence relation 


2F'(a;c; x) = aF(a + 1;¢; 2) — aF(a;c; 2), (11) 
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Eq. (10) may be written in the form 


(a. — a, — X) [ ea" "F(a, -c; 2) F(a. 3c; x) dz 


+ AA — 2a, +c — 1) [en *F(a, -c; 2) F(a, 3c; x) dx 


+ 2a. [ e*a°***F(a, 5c; 2)F(az + 130; 2) dz (12) 
| (a, — a, — A)F(a, 3c; x)F (az 5c; 2X) 
= ¢*y°**"" + a,F(a, 3c; x)F(a, + 1; ¢; x) 


— a,F(a, ;c; x)F(a, + 1; ¢; x) . 
In Eq. (12) we may let \ = 0, 1, 2, --- successively and thus obtain 
[ e *4°** "F(a, 5c; z)F (a, 3c; 2) dz for n= 60, 1,2, --- (13) 
provided that none of the resulting denominators is zero, which will be the case if 


a,— a, | #0,1,2, ---n. (14) 


Here also, a more symmetric recurrence relation for the integrals (13) may be obtained 
by interchanging v, and v, in Eq. (12) and subtracting the result from Eq. (12), giving: 


(a, — a) [ *0°**""F(a, 3c; 2)F (ae 3c; x) dx 


+ r(a, — ae) | e-*a°*** F(a, 5c; z)F(az 3c; x) dz 


+ Xa, / e*x°*** F(a, 5c; x)F(az, + 1;¢; x) dx 
(15) 
— da, | e*x°** F(a, 3c; x) F(a, + 1;¢; x) dx 
(a, — a.)F(a, 3c; x)F (az 5c; 2) 


=e?" + a.F(a, ;c; x)F(a, + 1;¢; 2) 





— a,F(a, 3c; 2)F(a, + 1;¢; 2)! 


from which the integrals (13) may be obtained by setting \ = 0, 1, 2, --- provided the 
condition given in (14) is satisfied. 

Two particular cases which may be obtained immediately from Eqs. (10) and (12) 
may be noted: 

1. If weset a, = 0,a, = ain Eq. (12), then, since F(0;c;z) = 1 (e ¥ 0, —1, —2, ---), 


ad 
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we have 


(a + 2) | e*x°*"'Fla;c; 2) dx — MA + — 1) [eco F Ga; ¢; x) dx (16) 
= —e72°""[(a + A)F(a;¢c; 2) — aF(a + 150; 2)] 

(a). Setting A = 0, 1, 2, --- successively in Eq. (16) we may obtain 
e*x°*""'F(a;c; x) dx for n= 0,1, 2, --- (17) 
provided that a # 0, —1, —2, --- —n. (It should be noted that for a = 0, 
—1, —2,--- —n, F(a; c¢; x) is a polynomial and hence the integrals (17) are elementary.) 
(b). On the other hand, setting \ = 1 — c, 2 — c, --- successively in Eq. (16) we 

obtain 

| e *x"F(a;c; x) dx for n=0,1,2,--- (18) 


provided c — a — 1 # 0, 1, 2, --+ n. 
(c). Finally, setting \ = —a, —a — 1, —a — 2, --- successively in Eq. (16) we 


obtain 
[ *2°~°-?-"F(a;c; 2) dx for n=0,1,2,--- (19) 
provided a ~ 0, —1, —2,--- , —(m — 1) andc —a — 1 # 0,1, 2, --- n. 


2. The other particular case to be noted may be obtained by setting a, = c,a, = @ 
in Eq. (10) and noting that F(c; c; x) = e*. We then obtain 


(c— a+) | a°'F(a;e; 2) dx + MA +e — 1) il a°*F(a;c; x) dx 


20) 
=2°"[(a+ + a)F(a;c; 2) — aF(a + 1;¢;2)]. 
(a). Setting X = 0, 1, 2, --- successively in Eq. (20) we may obtain 
[ x°*""'F(a;c; x) dx for n= 0, 1, 2, --- 21) 
provided a — c 0,1, 2, --- n. 
(b). Setting \ = 1 — c, 2 — ¢, --- successively we obtain 
[ x"F(a;c; x) dx for n=0,1,2,--:- (22) 
provided a ~ 1,2, ---n+1. 
(c). Finally, setting \ = a — c,a — c — 1, --+ successively we obtain 
x” *-"F(a;c; x) dx for n=z=Q0,1,2,--> (23) 


provided a —c + 1+ 1,2,--: nanda # 1,2,---n. 

IV. Bessel functions. As a final illustration we give a number of indefinite integrals 
involving the Bessel functions. Some of these do not appear in the literature; others 
that are to be found in the literature appear here as special cases. 
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Choosing Y, = J,(ax), a solution of 
2’Y" + 2Y’ + (a'r? — wp) Y = 0 (24) 
and y) = x*J,(8x) which satisfies 
ay!’ + (1 — 2d)ry’ + (62? + ° — v*)y = 0 (25) 


we have directly from Eq. (28) of [A], 
a J (ax) {2d2[2*J,(B2))’ + [@? — Br? +7 -— w? - d* JJ (Bx) } dx 
; (26) 
= x{J,(ax)[2J,(Bx)]’ — 2°J,(Bx) J i(ax)}. 
Here primes denote differentiation with respect to x, throughout. Using the recurrence 
relation 
aJi = —»J, + zJ,-; (27) 


Eq. (26) gives the recurrence formula 

[ 2-J,(ox {I(a? — B)x? + ( — d)? — p’]J,(Bxz) + 2dBrJ,_,(Bx)} dx 

x (28) 
= 2'[((\ —v + p)J,(ax)J,(Bx) + BrJ,(ax)J,-,(Bx) — atJ,-1(ax)J,(Bx)]. 

On the other hand, substitution of the recurrence relation 

gd. = pi) y= Bd vis (29) 

in Eq. (26) gives” 

| x (ax) { [(a? — B)x? + (v +r)’ - py |J,(B2) — 2rBxJ,+,:(Bx)} dx 

d (30) 
= 2A +» — p)J,(ax)J,(Bx) — Bard ,(ax)J,-,(Bx) + arJ,+:(ax)J,(82)]. 

Several interesting particular cases follow directly from Eqs. (28) and (30)’. 

1. Setting »y = B = 0, a = 1 in Eq. (28) we have 


| oy, dx + (dr? — yw’) | a" J,dz = al(a + uJ, — rJ,-1] (31) 
and setting \ = +u, +u + 2, tu + 4, --+ successively in Eq. (31) we may evaluate 
[a**J,de for n=0,1,2,---. (32) 

2. (a). Setting 
A—v= tu (33) 


in Eq. (28) gives 


r a 


(a’ — p’) | x" J (ax) J (Bx) dx + 2r8 | a J (ax) J,-,(Bx) dx 
: x (34) 


j 
/ 


= 2((\ —v + u)J,(ax)J,(Bx) + BaJ,(ax)J,-,(Bx) — axJd,-1(ax)J,(Bx)]. 





*Given in Ref. [1], p. 134, Eq. (7) for the particular case \ = 0. 
*Ref. [1], §5.1, p. 1383 Eq. (4). 
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If a? * 6’ and we let \ = 0, 1, 2, --- successively in Eq. (34) we obtain 
[ a*"'J,(aa)Ieys(B2) dz for n= 0,1,2, ++ (35) 


(b). If we let 
A+» = +u (36) 


in Eq. (30), we obtain 


(a? — 6°) / oy (ax) J,(Bx) dz — 2X8 | 2J,(ax)I».(B2) dx 


(37) 
= 2(\ +» — p)J,(ax)J(Bx) — Brd,(ax)J,+:(8x) + atJ,,.:(ax)J,(B2)]. 
If a? ~ 6’ and we let A = 0, 1, 2, --- successively in Eq. (37) we obtain 
[ a°'J.(e2)J.,-0(82) de for» = 0,1,2, ++ (38) 


The integrals (35) and (38) for a’ = #” may be found from the results for a’ # 6° by 
the use of |’Hospital’s rule. They are also expressible in terms of Bessel functions. 
3. Setting a = B = 1 in Eqs. (28) and (30) gives 


(v — »)* — a) [ PU, de + 20 il DJJ, de 


(39) 
=x —vt+pJid, + tJb,J 1 — tJy-adsl 
and 
(v +»)? — 7] [ 2"J,J, de — 20 / DI Sve, dt 
; (40) 
—_ x [(d + slices bd yd, = BS J 41 + tJ no18 ol 
respectively. 
(a). If we let \ = 0, —1, —2, --- successively in Eq. (39) or (40) we obtain 
[ e34. + & e0814%- (41) 
provided 
s+ een, 8-88 — 4) s+ mm (42) 


When one of the restrictions (42) is not satisfied, the integrals in (41) may be evaluated 
by applying |’Hospital’s rule to the results when the inequalities are satisfied and are 
expressed in terms of the functions 


’ 


| 2 
F(a) = S| Ga) = | (43) 
Ou - Ou iwwr 


‘Evaluated in Ref. [1], p. 135, Eq. (13) for the particular case n = 0, and for integer y, v and arbitrary 
integer n in Ref. [5]. 
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Actually, the integrals (41) will not depend on G, unless the factors multiplying the 
first integral in Eqs. (39) and (40) have a double zero, i.e. unless » = 0. If, however, 
un = 0 but » + O then the integrals (41) may still be expressed independently of G, , 
since they may be determined from the recurrence relations obtained by interchanging u 
and » in Eqs. (39) and (40). Thus, only if » = » = 0 will the integrals (41) depend on 
G, , and in this case they will depend only on G, . 

If we consider the integrals (41) for integer values of u, v, then we may note that 
differentiation with respect to v of the differential equation for J, : 


ey! + ad, + (2° — ¥)J, = 0 (44) 


gives 
°F!’ + 2F! + (x? — v’)F, = 2vJ,, (45) 
and therefore F,(x) may be expressed as a linear combination of the Bessel functions 
J,(z) and Y,(z).° Further, from Eqs. (27) and (29), we have upon differentiation with 
respect to v, 
2F) = oF, — zF,,, + J, (46) 


and 
xFi = —vF,+ 2F,_, — J, (47) 


respectively, and thus by induction F(x) may be expressed in terms of Bessel functions 
for any integer n. Thus for any integer values of y, v, the integrals (41) may be expressed 


in terms of Bessel functions, unless 1» = v = O (in which case G, is involved) and the 
-++)°, It may also be noted from all 


inequalities (42) are not satisfied (i.e. n = 0, 2, 4, 
this that for arbitrary u, v and regardless of whether or not (42) is satisfied, the integrals 
(41) may be expressed in terms of the solutions of the inhomogeneous differential equation 


ay” + ay’ + (2? — v)y = 2kJ,(z), (48) 


F being a solution for k = v, G, a solution for vy = 0, k = 1. 
(b). If weletAX =v tu,vtut+2,v+yu + 4, --- successively in Eq. (39) we 
obtain by recurrence 


oe dx for n=0,1,2,--- (49) 

provided 
yep XX —n, —n—1, —n — 2, +++ —2n. (50) 
Similarly, if we let} = —y tw, —yv t+u+ 2, —y+u + 4, --- successively in Eq. 


(40) we obtain 


a dx for n=0,1,2,--- (51) 


5In fact, Fo = wYo/2 (Ref. [1], p. 59, Eq. (5) and p. 64 Eq. (2)). 
6An evaluation for the particular case u» = v = 1, 2, 3, ... and n = 0 is given in Ref. [1], p. 137, 
Eq. (6) where it is noted that “no simple formula” seems to exist for the caseu» = » =n = 0. 
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provided 
—vtpHx —n, —n— 1, —n — 2, +--+ —2n. (52) 
With a slight change of notation the integrals (49) and (51) are expressible in the form 


| meee dx n= 0, Re 2, eee bu (53) 


It is worth recalling that in each of the integrals evaluated thus far we may substitute 
any solution of the differential equation being considered for the particular one that has 
been chosen here, provided only that the recurrence relation for this other solution is 
also substituted. Further, in connection with the integrals involving Bessel functions, 
since the differential equations and recurrence relations satisfied by J, and K, differ 
trivially from those satisfied by J, , we may also obtain integrals of the form given in 
(35), (38), (41) and (53) with J or K substituted for J. One particular case is worthy 
of note in view of the attention it has received in the literature®. Consider the integral 
similar to (53), with J, replaced by K, and J, replaced by either K, or J, . If we set 


u = 1/2 and note that, apart from a constant factor, K,,. = x '”e™*, then we obtain 
[arreK@ dx for n= 0,1, 2, --- (54) 
and 
[ 2° "e"L.@) dz for n=0,1,2,---, (55) 


the choice + in (53) giving the integrals (54) and (55) with n odd, —y giving those 
with n even. If, on the other hand, we had set a = —1, 8 = 1 in obtaining Eqs. (39) 
and (40) from Eqs. (28) and (30), then following the procedure used in deriving (54) 
and (55) we would have obtained 


[areK,@) dz for n = 0, 1, 2, --- (56) 


and 


[arve@) dx for n=0,1, 2, --- (57) 


Similar integrals, with the exponential replaced by sine or cosine and K and TJ replaced 
by J and Y may be obtained directly from (53) by setting u = 1/2. 

V. Extension to solutions of equations of arbitrary order. The procedure which has 
been used to evaluate in closed form integrals involving functions which are solutions 
of second order linear homogeneous differential equations may be extended to functions 
which are solutions of nth order equations quite simply. Equations (19) and (20) of 
[A] now become 


y Ee + b,(z) e | ou = 0 (58) 


i=0 





7Given in Ref. [1], p. 186, Eqs. (1) and (2) forn = 0. 
®The integrals (54) to (57) are given in Ref. [2] where it is noted that their closed form evaluations 


seems to have been overlooked in Ref. [4]. 
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and 


Y [a,(z) + 8b,(2)] = > b(2) 2 H(z, 0) Sn Ki, (59) 


i=0 
where Y is defined in Eq. (21) of [A]. Again using variation of parameters we may 
write in place of Eq. (4) of [A] 


fat T ¥e.. (60) 


where the Y, are linearly independent solutions of 


T's, (x) a = 0, (61) 
c(z) = a(x) + 8b,(z), (62) 
_, [ AK@ 
we / a, (63) 
| Y, Y, Y, 
| y! Y! y! 
| 
W = | 
(64) 
| 
} 
| 
| ae yor? a. ror 








= const. X exp E / = dx |, 


where 


-1) 


A; = cofactor of Y;""’’ in the determinant for W in Eq. (64), for n > 1, and (65) 
T9] 
A; = 1lforn = 1. 
Again, as with the second order equation, we assume that the boundary and initial 
conditions are such that y(z, t) = y(x, 0) = yo(x) so that 


- a (x) ~ 3 fe =0 (66) 


i=0 





and 
Y = 2 Yo (2). (67) 


The integrals 7; may now be determined from the n simultaneous equations which are 
obtained by substituting Eq. (67) in Eq. (60) and differentiating: 


LDYPL=y" k=0,1,2,---,n—1. (68) 
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In Eq. (68) no derivatives of J; appear in view of the equations 


> YY =0 k=0,1,2,---,n— 2, (69) 
i=l 
which may be verified directly from Eqs. (63), (64), (65). From Eq. (68) we may solve 
for the n integrals J; . We need consider only J, , however, since the solutions for the 
other integrals constitute merely a renumbering of the solutions Y; . The determinant 
for the J; is W, as comparison with Eq. (64) shows. We have, therefore, 


1,=W"D,, (70) 
where 
D, = the determinant obtained from W, Eq. (64), by replacing Y, by y . (71) 
From Eq. (7) of [A] and Eqs. (59) and (66) we have 


sK(z) = >> (c; — a,) d'Yo = a d'Yo (72) 


i=0 dx i=0 * dz’ 


and hence from Eqs. (63), (64), (70) and (72) we may write 


[ cx exp | / Cn=1 az |a, S¢.- a) Be | | Co-t az |p, (73) 
d Cy ‘=o dx Co 
which is the generalization of Eq. (28) of [A], or a similar equation with c; — a; replaced 
by ¢c; , which is the generalization of Eq. (30) of [A], and no longer requires that yo 
satisfy Eq. (66). 

It is seen that the idea developed in [A], combined with the extension presented in 
this paper, constitutes a simple and unified method for the evaluation in closed form of 
a large number of integrals involving the special functions. 
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—NOTES— 


ON ANALOGOUS THEOREMS OF FREDHOLM AND FRAME AND ON 
THE INVERSE OF A MATRIX* 


By DANIEL SHANKS (Naval Ordnance Laboratory) 


In this note we give a simple proof of Frame’s method of obtaining the characteristic 
equation of a matrix, and, by giving Fredholm’s method of solving his integral equation 
a simple non-determinantal formulation, we show that the Frame and Fredholm algo- 
rithms are identical. On the basis of this identity we investigate the possibility of 
approximating the inverse of a matrix by means of an abbreviated computation. 

Let A be an arbitrary N X N matrix and let 


AY + pry? + pry? + +--+ +p, = 0 (1) 


be its characteristic equation. Let 0 and J be the null and unit N X N matrices and 
define recursively the matrices, A, , and the numbers, a, , by: 


Ay = 0, Ao = & (2) 
A, = A(A,-; + a,-1]), a, = —* (trace A,). (n = 1, 2, +++). 


Then Frame shows (in [1], but in a different notation) that 

a, = Dr (3) 
and 

Aya, = 0. (4) 


A simpler proof may be adapted from the earlier work of Hotelling, Bingham and 
Girshick by using the theorem of Newton on the powers of the roots of a polynomial 


(2). If A, , Ax, «++ Aw are the roots of (1), and if 
8, = > rN, (5) 
P= 8, Po = —4(8 + Did) 
and, in general, n = 1, 2, --- , N, 
DP. = -| (8, + Di8a-1 + +++ + Da-18:). (6) 


Now for the matrix A, 
8, = trace A", 
and since, from (2), 
A, = A" + a,A"" + +++ +4,_,4, 


*Received February 8, 1954. 
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by computing the trace of both sides we see from (6) that Eq. (3) follows by induction. 
Further 
Ay + ayl = A” + a,A™" + +++ + ay1A + ayl = 0 
by the Cayley-Hamilton theorem. Then Eq. (4) is immediate. 
Since Frame shows that 


Ay + Gyl = A(Ay-1 + Gy-:I) + ay! = 0, 


the inverse of A is given by 
A™' = —(Ay-; + ay-,])/ay (7) 


if it exists. We now note that it also follows that if B = J — A then 


-1 7 A, 
Bt a1 4 204 (8) 


0 a, 


if it exists. 

Alternatively we may develop (8) as a corollary of the Fredholm theory of integral 
equations. This approach yields useful information concerning the rate of convergence 
of the two series in (8) and the possibility of approximating (8) by replacing the N 
there by some smaller integer, r. Fredholm’s solution of 


u(2) = fz) +f K(x, yuly) dy (9) 


may be written without determinants in a simple formulation which is virtually identical 
with (2) above. For from [3], Eqs. (15), (17) and (21) it may be seen that if we write 


K,(z, y) = 0, k, = 1, 


K,(z, y) = [ K(x, )K,-,(t, y) dt + k,,.K(z,y), (n= 1,2, ---) (10) 


= -" | K,(t, t) dt, 
N Je 


then the resolvent, 
o K,(x, y)" 


a: 


is a meromorphic function of \ such that Fredholm’s solution of (9) is 


R(z, y) = 


b 
u(x) = f(x) + / R(x, y, f(y) dy. (12) 
(It is perhaps worth noting in passing that if we replace k, = —(1/n)f? K,(t, é)dt in 
(10) by 
k, = 0 (n = 1, 2, ---) (10a) 


we obtain, instead, the “iterated” kernels K"”’ (x, y) and the Neumann series.) 
Now let a = 0, b = N, and assume that K(z, y) is a two-dimensional N X N step 
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function and that f(x) is a one-dimensional N step function. That is for 
i-ls2u<i 


(é,j = 1,2, --- ,N) 


j-lav<3 
we assume 
AK(a, y) = aj; (a constant) 
and (13) 
f(m) =f; (a constant). 


Then for A = (a,;) we find by comparison of (2) and (10) that K,(z, y)A” = A, and 
kd" = a, . Thus from (11) and (12), we see that for this degenerate kernel, K(z, y), 


Fredholm’s solution of (9) may be expressed matrix-wise: 


o- [ 4 Xe aby. (12a) 


0 a, 





However for this degenerate kernel we may also solve the linear system (9) directly: 
u = (I — A)"'f, (12b) 
and a comparison of the two solutions gives us (8) again. 


But this time we get more, for in proving that R(z, y, A) is a meromorphic function 
of \, that is, in proving the uniform convergence for all A, of the two series in (11), 


Fredholm shows that if 
| K(z, y) | S M*, 
then 
n/2 
| Ke, v) | $ py Mr - a. 

In matrix notation we thus find that if 

la; | S AM* = M, 
then the elements of A, are bounded in magnitude by b(A,), where : 

b(A,) = r'?M'N'/(r — DI. (14) 

Similarly 


la, | Sr°?M'N'/rl. (15) 


Therefore 





; "20MIN" (r + 1)°*”°(MN)""* ( MNe'” ) 
/n} 
| > G1 < > n”(MN)"/n! < CFD! mpg 


#=0 


If we use Stirling’s formula 


ee (<) 
al ~ (2rn)'”? \n 
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and let z = eMN(r + 1)7! this becomes 
| N 1 ntl 
| (og a 6 
2. od > [Qr(r + 1)]'*? 1 — ze” > 
Therefore if M is sufficiently small, that is, if z << 1 for some r < N, it becomes feasible 
to approximate (8) by 


=} 0 An 
B =I+ 2 A, . (8a) 


Do a, 
If we do so we are assured that the error in the denominator will be no larger than (16) 
and similarly the errors in the elements of the numerator will be no greater than 


(r+ abe Han 


N(2m)'? “1 = 2l(r + 1) /er]'” ° 
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THE ELASTIC SPHERE UNDER CONCENTRATED TORQUES' 
By ALFRED HUBER (University of Maryland) 


1. Introduction. In this note we consider an elastic sphere subject to two concen- 
trated torques acting at the end points of a diameter. The stress function occurring in 
this problem may be considered as a generalized axially symmetric potential’. This 
viewpoint enables us to obtain the solution by a simple potential-theoretic argument. 

The notation used in this article is the same as that used by Timoshenko in [6]. 

2. Stress function. We introduce cylindrical coordinates (r, 6, z) and consider the 
elastic sphere 2” + r’ < 1, subject to unit torques acting at the points r = 0,z = + 1. 
All relevant quantities appear in a meridian (z, r)-plane. According to the theory of 
torsion of circular shafts we must find a stress function (z, r) satisfying the following 
conditions: © is to be a solution of 

2 a2 
e425 to we (2.1) 


b 


Oz” or” r or 


defined and bounded in the half-circular region H = [r’ + 2” < 1] ( [r > 0], assuming 
the boundary value 0 on the open diameter, 1/27 on the open circular boundary. This 
problem shall in the sequel be referred to as problem (P). 

1Received May 1, 1954. This research was supported by the United States Air Force through the 


Office of Scientific Research of the Air Research and Development Command. 
For the definition of this concept the reader is referred to [8]. 
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We consider u(z, r) = 2r&(z, r) — 1. This function vanishes on the half-circle and 
can be continued analytically beyond it by the relation 


u(z, r) oe pu(zp’, rp), (2.2) 


where p = (z’ + r°)'”. For, by means of the correspondence principle (A. Weinstein 
[7]) the function v(z, r) = r~* u(z, r) is an axially symmetric potential in seven-dimensional 
space, considered in a meridian plane. The continuation (2.2) follows immediately from 
the classical reflection principle for harmonic functions, applied to v. 

u is now defined in the whole half-plane r > 0, satisfies (2.1), assumes known boundary 
values on r = 0 and u = O(p’) at infinity. Such a function is obtained as a generalized 
Poisson integral (A. Weinstein [8]) 


w= 34 f 
eee tL 


where d = [(z — ¢)’ + 7’]'”. An elementary computation yields 


“(eprtart [natar+ [ ratas}, ea 


— 5 [4 — (1 — g(a + b-*)\(2p' + 82 + 5r* + 2) 
(2.4) 
— 6e(a~* — b*)(1 — p')], 


where a = [(z + 1)? + r°']'”, b = [(e — 1)? + 7° ]'” and p = (2? + 7’). @ has the 
required properties. This is implied by the construction and can also be directly verified. 
We mention, that (2.3) could also be inferred from a formula of M. Olevskii [4], which 
was derived in a different way. Olevskii solved a more general problem whose solution is 
also accessible to the above potential-theoretic approach. Related results have recently 
been announced by P. Brousse [1, 2]. 
3. Stresses and displacements. The stresses are obtained by differentiation: 


eee or eS ee ee 
r or 4ar 
3zp” -3 -3 3 -5 -5 2\/o 4 2 2 
+ eo" — 6") + (a + b“)(1 — p)(2p° + 82" + Sr’ + 2) 
Tr Sarr 
+ ZU — pya* — &*) 
dar . 
a a _ 32 -3 -3\; 4 2 2 
ae 9s! ~eo [((a" + b“)\(p + 2p +2) 


+ (a~* — b™)(8zp" + 2)] — 3 (1 — p*)[((1 + 2a™* — (1 — 2b). 


Both components vanish on the axis of symmetry. In the equator plane (z = 0) the 
above expressions reduce to 


1. = : (2r° + 5° + 51+), 1 = 0. 
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Figure 1 shows the stress distribution in the equator plane, the dotted line indicating 
the respective stresses for the circumscribed cylinder. Figure 2 respresents the stresses 


in the cross-section z = — 0.5. 
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In the general case of the sphere z’ + r° < R’, subject to torques of magnitude M, 
the stresses at (z, r) can be computed by multiplying with MR~ the respective quantities 
for M = R = 1, taken at the point (zR™, rR™’). 

In order to determine the displacements for M = R = 1, we observe that u 
and v = r W(z, r) (see e.g. [6]), where V satisfies the equation 

ov , dv , 3av ' 
02 Or ' ror ° (3.1) 


=w= (0 
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and is connected with ® by the generalized Cauchy-Riemann equations 


ov O® 

Gr — = — 

" or’ 

G denoting the modulus of shear. V is determined up to an additive constant by these 
relations. Since Y = const. on z = 0 because of dV/dr = O6/dz = 0, we may define 


3;0¥ dF 
Gr tia (3.2) 


v(0,r) = 0. (3.3) 
It follows from (3.2) and (3.3) that 


_3 [f° 0® 
V(z,r) = Gr | Or dz. 
It can easily be seen, without carrying out explicitly the cumbersome integration, that 
v(z, r) can also be represented as a closed expression of elementary functions. 

4. Uniqueness of the solution. Let &, and ®, be two solutions of problem (P). The 
function f = ®, — ©, satisfies (2.1), is bounded in H and assumes the value 0 at all 
points of the boundary with the possible exception of A(—1, 0) and B(+1, 0). 

We first show that f approaches 0 also at A and B. Let f*(z, r) denote a solution of 
(2.1), defined in the quarter-circle Q@ = H (\ [z > 0], with boundary values f on z = 0 
(0 < r S 1), O elsewhere. f* exists; an explicit representation for f* in terms of the 
boundary values could easily be derived using the methods of Sec. 2. We now consider 
the function g = — f* which satisfies (2.1), is bounded in Q and vanishes on the 
boundary except possibly at B. By an inversion with respect to the circle of radius 1 
and center B, Q is mapped on Q, , an infinite subdomain of the half-plane r > 0. The 
function g, , defined in Q, by the relation 


g(z — 1,r) = b*g((z — 1b’, rb’), (4.1) 


where b = [(z — 1)? + r’]'”, satisfies (2.1) in Q, and vanishes on the finite boundary. 
Furthermore, g,; = O(p*) for p +. Applying the Phragmén-Lindeléf theorem in [3] 
to both g, and — g, , we conclude that g, = 0 throughout Q, . Hence g = 0 in Q. It 
follows that f approaches 0 at B. The same is true for A. 

Since f is continuous and identically 0 on the whole boundary of H we infer from the 
maximum-minimum principle that f = 0 in H. This proves the uniqueness. 

5. Discussion of singularities. Next we demonstrate that the solution of problem 
(P) yields the right type of singularity for the stresses at A and B. The state of stress 
corresponding to two concentrated unit torques in these points can be defined by the 
following limiting process. We subject the sphere to torsion by two homogeneous distri- 
butions of shearing stresses, of resulting moment 1 each, acting on the spherical caps. 
S, = [2 +7? = 1)N[|24+1| < 1/n) and S, = [2 + 7? = IO[|2—-—1]| < 1/n), 
n = 1, 2, 3, --- . Let ®, be the corresponding stress function, 73, and 77, the stresses. 
Their limits, as n ©, define the state of stress due to concentrated unit torques in 


A and B. 

Because of the maximum principle we have ®, < 1/2z, the boundary values being, 
continuous and satisfying the same inequality. By a similar argument we see that 
{@,(z, r)} is an increasing sequence. Therefore, it converges to a limit function &(z, r). 
The convergence is uniform in every compact subdomain of H and the derivatives of 


®, tend to the derivatives of &. It can be verified without difficulty, that @ is the stress. 
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function of the limit stresses and a solution of problem (P). Because of the uniqueness 
® is given by (2.4). 

This verification as well as the uniqueness proof are by no means superfluous. In 
fact, the solution ceases to be unique if the condition of boundedness is dropped. Con- 
sider for example the function 


@*(z,r)=r(1l—-2 —ry[e -—1)’? +r)”, (5.1) 


which satisfies (2.1) in H and assumes the boundary value 0 everywhere, except at B, 
where #* = O(b-’). This stress function corresponds to a singular self-equilibrated 
system of forces at B. Superposing (2.4) and (4.2) we obtain a function which fulfills all 
conditions of problem (P) except boundedness. In the terminology of E. Sternberg and 
F. Rosenthal [5], who treated the case of concentrated loads, it is a “pseudosolution”’ 
of the physical problem. 
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ON SUBHARMONIC SYNCHRONIZATION OF NEARLY-LINEAR SYSTEMS* 
By HIRSH COHEN (Carnegie Institute of Technology) 


When a linear system is driven by an external force, one has become accustomed to 
expect a frequency response related to the combination of the driving frequency and 
the natural frequency of the system. If one of the elements of the system is non-linear, 
but only slightly so, we may expect the frequency of the periodic response to behave 
like that of a linear system. For a physical system which is described by the non-linear 
differential equation 


i + f(y)y + wy = Fd, (1) 


where y is of the order unity, and ¢ is very much less than unity, this is true for all but 
@ small interval in the range of driving frequencies. It has been shown in several special 
cases, notably for the van der Pol equation [1], that if this system is driven by a force 
A cos (wt + ¢), the periodic response of the system for values of w? — w5 small, say of 
order ¢, is not a combination or beat frequency, but rather a function with frequency 









*Received May 10, 1954. 
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exactly that of the driving function, i.e. with circular frequency w. This is the phenomenon 
called synchronization (or entrainment, or locking-in). Various experiments involving 
synchronization are described in the treatise of Minorsky [2]. 

A further investigation [3] of the particular case, the van der Pol equation in which 
f(y) = — (1 — »’), has shown that a stable synchronized subharmonic response will 
also occur when the system is driven by a frequency near to three times the natural 
frequency of the system (“natural frequency” of the non-linear system refers to the 
natural frequency of the system for e = 0). Thus, when the driving force is of the form 
A cos (3 wt + ¢), and when w” — w, is small, a periodic response of circular frequency 
w will occur. 

It would seem natural to ask, then, under what conditions synchronization will appear 
for systems described by (1) for other driving forces. Suppose the system is driven by a 
function of the form A cos (nwt + ¢),— 


ij + ef(y)y + wy = A cos (nwt + ¢), (2) 


where A is of order unity and suppose, further, that f(y) is such that periodic solutions 
analytic in ¢ are known to exist (the conditions necessary for such solutions are given 
by Friedrichs who has used a method of Poincaré [4]). Under what conditions will sub- 
harmonic synchronization as described above be found? In the following discussion n 


will be considered greater than unity. 
The method to be used is a simple generalization of that given in [3]. The function 
y is expanded in a series in ¢, the first term being the “linearized”’ solution. 


y=Yyrtyuteyt-::: (3) 


yo will then contain two constants and these are determined to order ¢« by the casting 
out of secular solutions. A condition for the existence of synchronization and its interval 


can then be determined. 
Proceeding in this fashion, then, (2) is re-written in the form 


j + wy = A cos (nwt + ¢) + ef(y)y + &* — wo)y. 
Substituting the series (3), the equation of degree zero in ¢ is found to be 
tio + w'Yo = A cos (nwt + ¢) (4) 


and 
Yo = a coswt + A, cos (nwt + ¢), (5) 
where A/w (1 — n°’) = A, . The first degree in ¢ equation is 
ti + wy = —f(Yo)Ho + cyo , (6) 
where 
c = (w* — wi)/e. 


According to the method of Poincaré, it is now necessary to determine coefficients of 
cos wt and sin wt on the right hand side of (6), and to equate these to zero. The function 
f(Yo) is now periodic and may be expanded in a Fourier series, say 


F(Yo) = Ta, + - (A, cos mwt + b,, sin mut). (7) 


m= 
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Then for the right hand side of (6) we have 
cos wilca + 4wab, + 4wnA, sin o(a,-1 + Gns1) + 4wnA, cos o(b,., + B,-,)] 
+ sin wt[}wad, — Zwad, + fwnA, COS Y(Ay-1 — An41) — FunA, sin o(b,-1 — b,41)] 
+ terms in higher harmonics. 
Thus for c to be other than zero 
Swab, + ZwnA,[sin g(a,-1 + Qns1) + cos o(b,-1 + bn41)] 
must be other than zero. In order to obtain a and ¢ the coefficients of cos wi and sin wt 
are set equal to zero and solved simultaneously: 
ca + 4wab, + 4wnA, [sin g(a,_; + a,4:) + cos o(b,_, + b,.:)] = 0, (8) 
tay — aa, + 4nA,[cos g(a,-; — Gn+1) — sin g(b,-, — b,4:)] = 0. (9) 
This appears somewhat involved but in at least one important case, the coefficients 
Gz , bz , An—1 » On-1 » Gna1 » Onsi , are Such that these relationships become trivial. In the 
van der Pol equation, f(y) = —(1 — y). Expanding this with y. = @ cos wt + 
A, cos (n wt + g), the Fourier coefficients are found to be (for n > 3) 


ta, = —(1 — 40° — 44A)), a, = 4a’, b, = 0, 
a,-1 = aA, COs ¢, Qni1 = aA, COS ¢, 
(10) 
b,-1 = —aA,sing, b,-1 = —a@A, sin ¢. 
Inserting these quantities in (8) and (9) one obtains 
ca = 0, 4(1 — 4Aj) =’. (11) 


There is no range of synchronization, and “subharmonic resonance” (see [1] p. 323) 
occurs at w = w, only with amplitude a. In particular one observes that for A = 0, the 
known first approximation to the free vibration amplitude, a = 2, is assumed. 

For n = 3, the equations reduce to exactly those obtained in [3], but in this case, 
as was shown, there is a region of synchronization. 

It should be made quite clear that we have chosen a rather restricted definition of 
subharmonic response. A driving force of circular frequency nw might also produce a 
response of circular frequency say n w/r, where r is an integer less than n and this might 
also be termed subharmonic response. A further generalization can also be made in the 
damping term so that it takes on the form ef(y, y). The system is still of the nearly 
linear type and the analysis is quite the same. The equation is re-written 

ij + ef(y, 9) + 1r’w’y = A cos (nwt + ¢) + (r’w” — wi)y. (12) 

The series y = Y% + ey; + €Ye2 «+: is substituted in (12) and one obtains as conditions 
for existence of a synchronization interval 

c, = d,/a, , e, = 0, (13) 


where c,; = r’°w — w,/¢, a, = amplitude of response with circular frequency rw, and 


d, and e, are the Fourier coefficients of order r in the expansion of f(y, 7) as a function 
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of the zero degree solution, namely 
Yo = a, cosrwt + A, cos (nwt + ¢), A, = Ajw'(r’ — n’). 


It should be emphasized that the preceding work has been carried out on the assump- 
tion that the forcing term F(é) in (1) is of order unity. Synchronization effects have been 
sought for, therefore, only in the terms of the solution of this order of magnitude. 
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A NOTE ON THE PATTERSON FUNCTIONS* 
By C. R. PUTNAM (Purdue University) 


The following brief note was suggested by a paper of Hartman and Wintner [1], 
dealing with the Patterson function P(x), defined by 


P(2) = L" [ o(t)p(t + 2) dt, (1) 


where p(x) denotes a positive periodic function of period L. If p(x) possesses a finite 
number N of (relative) maxima on the interval 0 < x < L, then, as a consequence of 
what the above authors refer to as the Patterson principle in crystallography (see [2]), 
the function P(x) would have at most N(N — 1) maxima on the interval0 < x < L. 
(Obviously, the function P(x) is a positive even periodic function of period L, with 
absolute maxima occurring at the points x = 0, + L, + 2L, --- .) By constructing 
counter-examples for the case N = 2, the authors show that P(x) can have more than 
the predicted two peaks on 0 < x < L and conclude that the Patterson principle 
cannot be valid as a general mathematical theorem. 

The present note will deal with the case N = 1, so that p(x) corresponds, in the 
terminology of Patterson ([2], p. 521), to an electron density which is itself an “atomic 
function’’ possessing a single peak on 0 S x < L. It will be shown that not only is the 
Patterson principle (which would here deny the existence of any peaks of P(x) on 
0 < x < L; see loc. cit., middle of page 522) false even in this case, but that P(x) may, 
in fact, have any specified finite number of peaks on 0 < x < L. 

To this end, let n denote an arbitrary positive integer and divide the interval 
0 < x S L into 2n parts of length d = L/2n. Next define the step-function p(x) 


*Received May 24, 1954. 
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{= p(x + L)) on0 S 2x < Lby: 
1, Os 2 < d, 


II 
~) 


p(x) = 9k, (2 — 3)d S x < (2k — ljd (k = 2, +++ ,n), (2) 


Lb, (2n — l)d S x < 2nd, 


where b is a positive constant to be determined. Then the graph of P(x) is a polygonal 
path with vertices having abscissas on 0 S x < Latx = 0,d,---, (2n — 1)d. It will 
be shown that for a proper choice of the constant b the function P(x) of (1) has maxima 
occurring at the points x = 0, 2d, --- , (2n — 2)d. For reasons of continuity, it is clear 
that the above p(x) can be replaced by a continuous positive function possessing a single 
maximum on the interval0 < x S L and so that the corresponding P(x) of (1) possesses 
n maxima on 0 S x < L; see the remark of [1], p. 272. 
An easy calculation shows that, for the p(x) of (2), the function P(z) satisfies 


| (n + 1)bd, for m=1,3,°---,2n—l, 
P(md) = c, + 4 (3) 
\(n + 2)bd, for m = 2,4, --- ,2n — 2. 
Accordingly, since the constants c,, do not depend on ), it is clear from (3) that for a 
sufficiently large b, the inequalities P(2kd) > P[(2k + 1)d] and P[(2k + 1)d] < P[(2k + 
2)d] hold for k = 0, 1, 2, --- , nm — 1. Thus, the function P(x) has n maxima 
on 0 S zx < L as claimed, and the existence of a continuous p(x) with the desired 
properties now follows as remarked above. 
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ON HELICAL SPRINGS OF FINITE THICKNESS* 
By PETER HENRICI (National Bureau of Standards) 


1. Introduction. The problem of determining the stress distribution in a statically 
loaded, closely coiled helical spring, whose cross-section is not necessarily small in com- 
parison with the diameter of the helix, has been solved approximately by iteration by 
Gohner [5, 6] and exactly in terms of a series of appropriate Legendre functions by 
Freiberger [4]. In the present paper it is shown that the exact values of the stress con- 
centration factor as well as of some other significant quantities of the problem can be 
developed in terms of a certain parameter (namely the reciprocal of the so-called spring 
index) so as to recover Géhner’s approximations and to obtain further approximations 
of the same type. 


*Received May 20, 1954. The preparation of this paper was sponsored by the Office of Naval 


Research. 
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2. Formulation of the problem. Mathematically the stress problem under con- 
sideration may be stated as follows’. Let the z-axis of a system of cylindrical coordinates 
(r, 8, 2) coincide with the spring axis and denote by G the intersection of a turn of the 
spring with the half-plane # = 0. Let the center of gravity of S lie on the plane z = 0 at 
a distance R from the z-axis. Let furthermore 

G sin 2a 


Co am, 


where G is the shear modulus of elasticity and a the difference of the pitch angles of the 
spring in loaded and unloaded state’. Then, if the “stress function” ¢(r, z) solves the 
boundary value problem 





a°o 
or” 





2 
oF 32 1g 


+= 
Oz r or (1) 


¢=0 __ onthe boundary of SG, 


the shear stresses (7, , 7,) in G are given by 


2 2 
52 o%,&. (2) 


fu — 
r° dz’ 4 


r= —C 


Here C, because it involves the unknown quantity a, is unknown. If one wishes to 
compute the maximum shear stress occurring under a given load, 


Op 2 ag 2 41/2 
eee (nae 
max C Tmax b J Tmax M AX (R r) | or + Oz ’ (3) 
one has to eliminate C from (3) and from the expression for the total moment,’ 
W = CW’, 
, 2 Fil dp dd 3 I -3 
* — —r—-2z>— g= —2 
W R || r | ce 7 “s i dr di R r “ dr dz, (4) 
Ss S 


for which by mechanical considerations one can also write 

W = RF, 
where F denotes the force exerted on the spring in the direction of the z-axis. It follows 
that 





a 7 
won > RP | | 
T R we | 

3. The stress-concentration factor. We are here primarily interested in the case 
where G is a circle. If we denote its radius by p, it is convenient to write 


Tmax = K. rn ’ (6) 


(5) 


and Goodier [10], pp. 391-395; Wahl [11], pp. 38-42; see also Collatz [3], pp. 292, 310, 313, 355, 363, 
where the problem is used as an example for the application of various numerical methods. 

The theory is exact only if the pitch angle of the spring is zero, i.e., if the spring reduces to an 
incomplete tore. However, as the quoted authors point out, it may be regarded as a good approximation 
in the case of a small pitch angle. For a detailed discussion of this point see Wahl [11], p. 41. 

’Mitchell [8], p. 134. 
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where 
0) 2RF - 
Tmax = ran wy (7) 

Tp 
is the maximum shear stress as given by the well-known elementary theory (neglecting 
the finite thickness or, what amounts to the same, the curvature of the spring), and 
K, is the so-called stress-concentration factor*. Since in this case it is clear for mechanical 
reasons that r,,,, occurs at the point P of the circumference of the circle nearest to the 
spring axis, we obtain from (5), (6) and (7), denoting the shear stress at this point by 
7 i= Crt ’ 
rp | tt 
, P / 

K, = “9 | we!’ (8) 
4. Solution of the boundary value problem for a circle. As indicated by Mitchell 
[8], an explicit solution of the boundary value problem (1) can be given for the case 


that G is a circle if bipolar coordinates (£, 7) defined by 


c sinh é ssinh 7 
eS oT 


cosh § — cos 7’ “cosh £ — cos 7 


(c = const.) are introduced. It is easily seen that if & > 0 is a fixed value of ¢ and 7 
varies from 0 to 27, the point (r, z) travels along a given circle with center at (0, R) and 
radius p < RF provided 

cosh & = R/p, c = (R’ — p)"”. (9) 
For the following it is convenient to denote by 

q = p/R 

the recipro¢al value of the so-called spring index. The solution of (1), which in essence 
was first given explicitly by Freiberger [4] (who, however, uses a stress function slightly 
different from ¢), can then be written as follows: 


g(r, 2) = 2°°'r"'(R’ — p’)(cosh é — cos n) *”’ sinh’ & 


. >> €,{ A,(q)Q;-1/2(cosh £) — Q,-1/2(cosh £)} cos nn. 
Here e, is the Neumann number (¢, = 1; ¢, = 2, n = 1, 2, 3, ---), 
A,(g) = Qn-ie(q )/Qr-rA(q); (11) 


and Q,(z) and Q*(z) denote the ordinary and the associated Legendre function of the 
second kind as defined in [7]. It can be shown that the series in (10) converges uniformly 


in £ and 7 for all 7 and all ¢ such that cosh""§ S gq < 1. 
From (2) and (10) one readily derives the following expansion for r#, which is equiva- 


lent to that given in [4]: 


rt _ ry ae ie q) ‘(1 + q)'”” 
° (12) 


Do (—D"ef(m — 8/2) An(DQ+12(97') — (a + 1/2)Qns12(9"')}. 


n=0 


‘Wahl [11], p. 110. 
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Similarly one finds from (4) and (10) 
Wwe =we?s+ we”, (13) 
where 
We = weg “[l — (1 — 9°)", (14) 


whilst W*‘ is given by the more complicated expression 


We = 9 pg (1 — — (15) 
S.(q) = a (s? — 1)""Qy-1/2(8)Q2-1/2(8) ds. (16) 


These formulas can again be shown to be equivalent to those given by Freiberger. 
Formula (15) has the numerical advantage over the corresponding one in [4] of not 
involving differences of large numbers’. 

5. Developments in powers of g. In order to bring the above results into a practical 
form, we develop the right sides of (12) and (13) in terms of powers of the only essential 
parameter g. This is done by developing the Legendre series termwise and rearranging. 
The rearrangement, which is permissible in view of Weierstrass’ theorem, is facilitated 
by the fact that the nth terms in (12) and (15) are 0(g") and O(q*"*”) respectively, so 
that in collecting powers of a fixed order one has to cope only with a finite number of 
terms. Computational details are omitted. 

We start by developing A,(qg). Using (11) we obtain by straightforward division of 
the two hypergeometric series involved in the definition of the Legendre functions 








Py. 1g t-ch 2 _12n'+8&m-3 
aNd = a + 1/2)(n + 3/2) n+1_% ~ 16m + In + 2) 2 
(17) 
_ 80n* + 288n* + 8n* — 392n + 33 5 _ sob 
128(n + 1)°(n + 2)(n + 3) 
It is now easily found that 
a=? 215 5 , 6827 sb, 
ac elit eat ia + 3a 7 + gig2 2 + (18) 
For W*'” the coefficients of the power series can be given in closed form: 
/2 ‘. 
we in SBS en (1/2), _ 2). ” ag (19) 


(n+ 1!? 


In order to expand W*’, we need to develop the integrals S,(q¢). By partial integration 
one finds 


8) = FS QanOP |)  — f=) (aaOT a. 


1 


- 


‘See the remark in [4], p. 370. 
8(a)o = 1; (an = a(at+1)...(€a+n-—1),n =1,2,.. 
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The terms of this expression are easily expanded by use of the formula 


" tea a _— 1)~*[Q3(s)]? 


~ pinef TG a + pra/ay + ut Dem p Pe agptinch’, ni 20) 
= € o*'Te + 3/2) pan n! 34 2 Dy f. 3, » rs 3/2 . 


which can be proved by application of a result due to Clausen’. Carrying out the ex- 
pansions and rearranging, we finally find 


4 
reco oP J, 1 ie 1 _ 849 +} 2 
W 2¢ { 4% +397 — d4576 2 + a 








and hence, using (13), 
a 
rs _ me f 2269 
W 4 \! + 


ts 
1 + iopgg ¢ + vf (22) 


16 








In virtue of (8) we obtain now for the stress-concentration factor the expression 


eve tral eatin sMibes.:. " 
Ke=1+qqd+ 9 +9569 + o45767 + °° (23) 





6. Comparison with other approximate results. The first three terms of the last 
formula, as well as the corresponding terms in (18) and (22), agree with Géhner’s second 
order approximation [5, 6], which was constructed by carrying out the first two steps of 
an iterative method of solving (1). See also [2] and [10] for expositions of the method. 
The third order term in (23) however does not agree with a term + q* given by Wahl’. 
The error seems to be due to the fact that, following [10], Wahl writes (in our notation) 

pj 1 1 1 »\ 
ton * git theta}, 


~ 


(which up to the term of second order in gq is identical with (18)) and then extends the 


expansion only of (1 — q)~’ to the third order term. 
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BOOK REVIEWS 


Wave motion and vibration theory. Edited by A. E. Heins. Proceedings of Symposia in 
Applied Mathematics. Volume 5. McGraw-Hill Book Co., Inc., New York, 1954. 


v + 169 pp. $7.00. 


The volume contains the papers presented at the Fifth Symposium in Applied Mathematics of the 
American Mathematical Society held at the Carnegie Institute of Technology in June 1952. The four 
sessions of the Symposium were devoted to Stability of Fluid Motion, Hydrodynamic Waves, Diffraction 
and Scattering Problems, and Vibration Theory. The contributors to the volume are: C. C. Lin, 8. Chan- 
drasekhar, P. R. Garabedian, W. Bleakney, N. W. McLachlan, H. Levine, W. Magnus, A. E. Heins 
and H. Feshbach, G. F. Carrier and W. H. Munk, J. J. Stoker, E. W. Montroll and J. M. Greenberg, 


E. H. Lee, A. Weinstein, 8. Lefschetz, and R. J. Duffin and A. Schild. 
W. Pracer 


Relaxation methods. By D. N. de G. Allen. McGraw-Hill Book Company, Inc., New 
York, 1954. ix + 257 pp. $7.50. 


The aim of the present volume is clearly and accurately stated in the author’s preface: “... to 
explain as clearly as possible how to relax; its contents were entirely chosen with this end in view... . 
not intended to be a reference book ... no attempt to discuss the more advanced and difficult applica- 
tions [such as] non-linear equations. Nor . . . to discuss the physical facts and assumptions. . . Another 
volume could be written dealing with such topics as convergence, finite difference and truncation errors, 
parallel procedures used in relaxation and other computational methods, and so on... . I have deliberately 
kept off these subjects. . .” 

When an author defines his aim, the duty of the reviewer is to discuss the extent to which he has 
achieved it; the prospective reader or purchaser can make his own decision as to whether his purposes 
will be satisfied by the author’s aim. Viewed in this manner, the task of reviewing the present volume 
is a pleasant one. Mr. Allen writes with great clarity. The book is well arranged, progressing from simple 
to difficult problems. The chapters are reasonably self contained so that a reader already familiar with 
the beginning subject matter may turn at once to the later chapters. Although the author does not 
discuss in detail any physical problems as such, he does not hesitate to draw upon physical analogies 
to explain portions of the relaxational technique. An important example of this is in connection with 
Poisson’s equation with curved boundaries where the relaxational net is viewed as a physical entity, 
being an approximation to a continuous membrane. 

The detail of the book appears to be ideal for the serious student. Beginning examples are worked 
out in full, and the new parts of each new concept are fully illustrated. Once the problem has been reduced 
to the point where further detail would be repetitious the remainder is left as an exercise. Frequently 
the final result is quoted so the reader can check his result. Finally, each chapter concludes with a few 
problems. As the author points out in his preface, “The examples at the end of each chapter have been 
designed to give the student . . . practice; if they seem to be uninspiring drill problems, well, so they are. 
No relaxer can learn too early to take the rough with the smooth. . .” 

The first four chapters of the book are devoted to linear algebraic equations, frameworks, and 
ordinary differential equations. Next come six chapters on partial differential equations including 
Laplaces’s, Poisson’s, quasi-plane-potential, biharmonic, and simultaneous. Also discussed in these 
chapters are curved boundaries, normal-gradient boundary conditions, and triangular nets. Next come 
two chapters on eigenvalue problems and two on internal and free boundaries. The penultimate chapter 
gives an ingenious technique by which the relaxation method can be used to solve non-elliptic differential 
equations. Finally, the last chapter applies the relaxation method to an example of Laplace’s equation 
in three dimensions. The book then closes with four appendices, two of which summarize useful formulae 
and two of which give mathematical details omitted in the text. 

Despite the blurb on the jacket, this is the second book on relaxation, the first being “‘An Introduc- 
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tion to Relaxation Methods” by F. S. Shaw (Dover, 1953). Both books are written by former students 
of Sir Richard Southwell, and both are approximately the same length and cover approximately the 
same ground. Shaw’s book is entirely mathematical, while as noted Allen’s makes some use of physical 
analogies. Shaw treats complex examples in considerably more detail, but does not contain any problems 
for the reader. On the other hand, Shaw has nothing corresponding to Allen’s last two chapers. If a 
conclusion is necessary, this reviewer would probably prefer to teach a course using Allen’s book as 


Xt, bu woul e€ a copy or eacn h 8 Own DOOKCase. 
t t ] ] ] k f I in 1 ) I v k ise 
P. G. Hopce, JR. 


Four-place tables of transcendental functions. By W. Fligge. McGraw-Hill Book Co., 
Inc., New York, and Pergamon Press Ltd., London, 1954. 136 pp. $5.00. 


According to a statement in the preface, these tables were prepared for ‘“‘those who are old-fashioned 
enough to do their daily chores with the slide rule, but modern enough to use mathematics beyond sines 
and cosines whenever that will help solve a problem.” The useful volume contains tables, of slightly 
more than slide rule accuracy, of those transcendental functions that are most likely to occur in problems 
of physics and engineering (trigonometric and hyperbolic functions; exponential function and logarithm; 
Bessel functions; Thomson functions; elliptic integrals; error function; Fresnel integrals; sine, cosine, 
and exponential integral; gamma function). The volume also contains a collection of formulas concerning 


these functions. 
W. PRAGER 


Stochastic processes. By J. L. Doob. John Wiley & Sons, Inc., New York, and Chapman 

& Hall, Ltd., London, 1953. vii + 654 pp. $10.00. 

This book contains the first complete and detailed treatment of stochastic processes and is essential 
for anyone concerned with the mathematical] basis of the subject. 

The first two chapters deal with the necessary probability background and definition of the various 
classes of stochastic processes—Gaussian, Markov, stationary, martingales, etc. There follow detailed 
sections on processes with mutually independent, uncorrelated and orthogonal random variables, discrete 
and continuous parameter Markov processes, martingales, processes with independent and orthogonal 
increments and discrete and continuous parameter stationary processes. The final chapter treats the 
theory of linear least squares prediction. 

There is a supplement at the end of the book in which the author outlines those aspects of measure 
theory most needed for the subject matter of the book. However, a fairly thorough knowledge of measure 


theory may nevertheless be considered to be an essential prerequisite for the profitable use of this work. 
LEonaRD C. MaxIMon 


Theory of matrices. By Sam Perlis. Addison-Wesley Press, Inc., Cambridge 42, Mass., 

1952. xiv + 237 pp. $5.50. 

The book gives a concise and clear treatment of the theory of matrices with emphasis on the basic 
ideas rather than particular applications. It can serve well as a text for a year course, especially if supple- 
mented by applications of the material presented in the book. The approach of the text emphasizes 
strongly the algebra of matrices, with chapters on vector spaces, equivalence, rank and inverses, congru- 
ences and Hermitian congruences. Other topics considered include determinants, polynomials over a 
field and matrices with polynomial elements, similarity, characteristic roots and vectors, and linear 
transformations. Exercises are given at the end of each chapter. 

The book provides a good background for anyone who is to specialize in either pure mathematics or 
physics and applied mathematics. Although some prior acquaintance with the basic notions of modern 
algebra would help the student using this text, the concepts are not presented in a form so abstract that 


the student without such background will feel lost. 
LEONARD C. MAxIMON 


























